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ABSTRACT 

We  will  show  that  utilizing  precursors  significantly  improves  the  preformance  of  an  HRR  radar  system  for  foliage 
penetration.  We  will  begin  by  explaining  the  many  advantages  of  HRR,  and  its  free-space  robustness  to  noise.  We 
will  explain  the  limitations  of  this  standard  system  in  a  highly  absorbing  and  dispersive  media  such  as  foliage.  We  will 
conclude  by  outlining  the  design  of  a  precursor  based  ERR  system  which  significantly  outperforms  a  standard  HRR 
system  in  the  presence  of  absorption  and  dispersion.  This  new  system  seems  to  be  viable  for  foliage  penetration. 

!•  Introduction 

Identifying  targets  under  a  foliage  canopy  has  proven  to  be  a  very  difficult  problem  for  the  military.  The  current 
trend  for  many  radar  systems  is  the  design  of  high  frequency  systems  with  high  resolution.  These  systems  do  not 
perform  in  the  presence  of  a  highly  scattering  and  absorbing  media  however.  High  Range  Resolution  (  HRR  )  radar 
systems  have  many  positive  characteristics,  however.  They  utilize  multiple  pulses  to  significantly  improve  upon  the 
resultant  signal  to  noise  ratio  (  SNR  )  of  the  system.  In  addition,  the  mixers  in  these  systems  produce  constant 
phase  outputs  for  each  pulse  returned  from  a  theoretic  point  target.  This  constant  phase  pulse  is  then  integrated 
over  the  length  of  the  pulse,  additionally  increasing  the  final  SNR  of  the  system.  The  basic  design  of  an  HRR  system 
does  not  account  for  absorption  and  dispersion.  We  will  show  that  the  performance  of  such  a  system  is  significantly 
degraded  in  a  foliage  type  setting. 

The  goal  of  this  paper  is  to  explain  how  a  hybrid  system  should  be  designed  to  utilize  the  advantages  of  HRR, 
while  minimizing  the  negative  effects  of  dispersion  and  absorption.  This  system  is  fundamentally  designed  around 
precursors. 

Precursors  were  first  explained  in  1914  by  Sommerfeld  and  Brillouin.^’^  Their  primary  focus  was  on  analyzing 
these  pulses  with  respect  to  causality.  The  conclusion  of  these  papers  was  that  the  leading  edge  of  the  Sommerfeld 
precursor  travels  at  the  vacuum  spe^  of  .light,  but  is  absolutely  causal.  Similarly,  the  Brillouin  precursor  tends  to 
have  a  majority  of  its  energy  at  the  leading,  and  corresponding  trailing  edges  of  the  pulse  and  is  causal.  Unfortimately 
Brillouin  concluded,  due  to  incorrect  asymptotics,  that  these  pulses  would  not  be  useful.  Recent  work  has  shown 
that  Brillouin’s  analysis  was  flawed  and  that  the  Brillouin  precursor  is  the  dominant  term  in  the  far  field  of  any 
dispersive  and  absorbing  media.  ^ 

The  fundamental  property  which  makes  precursors  so  attractive  is  that  they  are  not  exponentially  attenuated  by 
the  media  through  which  they  pass.  The  energy  in  a  precursor  is  only  attentuated  at  a  rate  of  1/ y/z,  as  opposed 
to  €xp{—ct{fk)z)  for  a  given  frequency  fk  and  absorption  coefficient  Oi{fk)*  We  will  show  that  the  utilization  of  this 
property  is  fimdamentally  important  if  one  is  to  acquire  significant  signal  strength,  through  an  absorbing  media 
layer. 

This  work  was  supported  by  the  Air  Force  Office  of  Scientific  Researches  Electromagnetics  Program. 
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2.  Standard  HRR  Processing 

W^hen  trying  to  design  new  radar  systems,  it  is  very  important  to  understand  the  positive  aspects  of  the  gold  standards 
which  are  utilized  in  similar  applications.  This  is  especially  true  when  the  goal  of  the  new  system  is  far  more  difficult 
than  that  of  the  standard  systems.  Synthetic  High  Range  Resolution  Radar  (  HRR  )  has  not  been  succesrful  in  a 
foliage  penetration  situation,  but  it  is  the  standard  among  one-dimensional  processing  algorithms.  In  addition,  the 
gold  standard  in  two-dimensional  processs,  Synthetic  Ai>erture  Radar  (  SAR  ),  utilizes  many  of  the  same  techmques. 
We  will  briefly  outline  the  positive  characteristics  of  HRR,  and  the  shortcomings  of  HRR  so  that  we  can  focus  the 
direction  of  this  program. 

HRR  radar  systems  are  commonly  used  for  both  attack  and  surveillance  by  the  military.  Radar  systems  engin^r- 
ing  is  often  times  a  tradeoff  between  achieving  proper  resolution  and  maintaining  acceptable  signal  to  noise  ratios. 
Synthetic  HRR  is  a  technique  to  overcome  these  difficulties. 

The  basic  idea  is  to  transmit  N  different  signals,  which  are  N  times  longer  than  the  desired  resolution,  and  use 
the  N  different  responses  to  recover  the  desired  resolution.  There  are  several  advantages  to  this  approach,  from  a 
signal  to  noise  point  of  view. 

To  begin  with,  we  will  assume  that  the  HRR  system  is  transmitting  at  N  different,  equally  spaced  frequencies,  /*. 
Thus  the  polarized  signal  can  be  represented  as  yk{i)  ~  (cos(27r/fci)  ism{2irfkt))  X^{t)  —  exp(27r«/fct)A'^(t),  where 
X^(t)  is  a  characteristic  function  representing  the  time  length  of  the  signal.  The  received  signal  will  arrive  with  a 
time  lag  r  =  2R/c,  where  R  is  the  distance  to  the  target,  and  c  is  the  speed  of  light.  Thus  neglecting  target  velocity, 
the  returned  signal  can  be  represented  as  Vkit)  =  A^xp{27vifk{t  -  r))X\t  -  r),  where  A  is  a  constant  proportional 
to  the  RCS  of  the  target. 

Upon  reception  this  signal  is  mixed  with  the  matched  filter,  Zk{t)  =  exp(— 27ri/fct),  and  the  resulting  output^ has 
constant  phase  and  amplitude,  and  thus  the  output  of  the  mixer  is  a  constant  phase  signal  Ok{t)  =  exp(“27ri/fcT)A'  {t- 
t).  At  this  point  in  time,  it  would  be  p)ossible  to  solve  for  t,  by  standard  arithmetic  means,  such  as  an  inverse  tangent. 
It  could  also  obviously  be  measured  by  the  X\t  -  r)  term.  There  are  two  signal  to  noise  improvements  which  make 
HRR  more  signal  to  noise  robust  than  these  basic  techniques  would  impart. 

The  first  is  that  by  making  the  signal  N  times  longer  than  the  desired  bandwidth,  one  can  easily  increase  the 
signal  to  noise  ratio  by  integrating  or  averaging  the  constant  phase  signal  Ojfc(t),  which  will  increase  the  signal  power 
by  a  factor  of  iNT,  while  only  increasing  the  noise  power  by  a  factor  of  y/N,  Thus  the  signal  to  noise  ratio  had  been 
increased  by  a  factor  of  y/N.  Once  again,  at  this  point,  you  could  directly  estimate  t ,  as  long  as  the  possible  range 
of  r’s  respects  the  periodic  nature  of  phase. 

HRR  takes  this  process  one  step  further,  however.  The  system  transmits  N  of  these  signals  g/fc,  and  measures  the 
averaged  phase  values  of  each.  Thus  N  values  Ofe(t)  =  exp(-27ri/fcT)  are  measured  and  stored  in  a  vector  6.  Now, 
assuming  that  the  frequencies  of  the  fk  axe  evenly  spaced,  and  neglecting  discretization  effects  of  the  continuous 
variable  r,  this  vector  6  can  easily  be  seen  as  a  frequency  shifted  column  of  a  standard  FFT  matrix.  Thus,  an  inverse 
FFT  will  result  in  an  approximate  5  function  response,  at  the  time  lagged  position  r . 

The  second  signal  to  noise  effect  is  achieved  above.  The  orthogonal  matrix  takes  a  constant  amplitude  vector, 
and  transforms  all  of  the  power  of  that  vector  into  one  fine  range  bin,  increasing  the  localized  signal  power  by  a 
factor  of  N,  Assuming  that  the  noise  is  white,  however,  the  noise  power  will  be  evenly  spread  across  the  range  bins, 
so  the  signal  to  noise  ratio  has  been  increased  by  a  factor  of  N  in  this  step.  The  total  increase  in  signal  to  noise,  in 
this  idealistic  model,  is  therefore  Ny/N.  This  process  is  illustrated  in  Figure  1. 

2,1.  Absorption  Problems  and  Solutions: 

The  performance  of  HRR  in  freespace  is  well  documented.  We  are  interested  in  penetrating  an  absorbing  and 
dispersive  media,  however.  A  number  of  the  assumptions  in  the  above  outlined  HRR  discussion  break  down  in  this 
case.  First,  the  returned  signals  will  not  have  constant  phase  when  they  exit  the  mixer.  More  importantly,  they  will 
not  have  equal  amplitude,  but  rather  will  be  absorbed  at  different  rates.  Finally,  much  less  power  will  return  to  the 
receiver,  due  to  the  absorbtion.  Much  of  the  noise  in  the  system  is  signal  independent  (  ambient  noise  ),  so  the  SNR 
will  be  greatly  degraded. 

We  will  now  outline  the  basis  of  our  precursor  based  HRR  system.  This  system  is  explained  most  easily  with 
the  help  of  Figure  2.  In  order  to  induce  the  precursor  phenomena  we  swap  the  phase  of  the  real  and  imaginary 
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IMoisse  level  is  S  dB 


Figure  1:  The  robustness  of  a  standard  free-space  HRR  system  is  illustrated  above.  A  point  target  was  theoretically 
placed  at  the  point  of  the  first  peak.  The  system  transmitted  128  pulses  that  were  128  feet  long.  The  processing 
outlined  above  produced  a  final  resolution  of  1  foot.  The  noise  level  at  the  receiver  was  5  dB,  and  the  noise  is  almost 
completely  suppressed.  The  secondary  peak  is  a  result  of  ’’range  walk”,  to  the  second  gross  range  bin.  For  foliage 
penetration  this  effect  is  negligible,  because  the  ground  will  act  as  a  barrier,  and  objects  outside  of  the  128  foot 
initial  gross  range  bin  will  be  ignored.  Without  dispersion  or  absorption,  this  system  will  perform  well  even  at  noise 
levels  of  -10  dB. 

components  of  the  transmitted  signals,  at  every  other  zero  crossing.  An  additional  detail,  however,  is  that  the 
real  and  imaginary  components  are  phase  lagged,  from  a  standard  exponential.  Note  that  the  real  and  imaginary 
components,  shown  in  the  top  of  Figure  2,  are  tt  out  of  phase  with  each  other.  This  is  necessary  to  make  sure  that 
the  received  signals,  shown  at  the  bottom  of  Figure  2,  are  7r/2  out  of  phase,  as  would  be  the  case  with  a  standard 
exponential. 

Precursor  HRR.:  The  basic  steps  in  the  precursor  HRR  system  are  outlined  below.  Note  that  with  the  exception 
of  the  transmitter,  everything  else  is  nearly  identical  to  traditional  HRR. 

•  Transmit  signals  that  are  phase  swapped  every  cycle,  and  phase  delayed  as  illustrated  in  the  top  illustratiion 
of  Figure  2. 

•  The  signals  travel  through  the  material,  and  are  received  as  in  the  bottom  illustration  of  Figure  2. 

•  Mix  the  received  signals  with  appropriate  exponentials  producing  nearly  constant  phase  outputs  as  in  HRR. 

•  Integrate  the  output  of  the  mixer  over  appropriate  large  range  bins. 

•  Collect  128  of  these  signals,  and  use  the  matched  filter  to  produce  the  final  output. 

To  give  a  baseline  for  our  material  model,  we  illustrate  its  effect  on  one  of  the  standard  HRR  signals.  In  Figure 
3,  we  have  transmitted  a  standard  HRR  signal,  and  its  corresponding  precursor  HRR  signal.  Both  signals  originally 
were  at  the  same  frequency,  with  phase  swapping  used  on  the  precursor  HRR  signal.  Note  that  the  precursor  signal 
was  received  with  5.8  times  more  power  than  the  standard  HRR  signal.  This  effect  gets  larger  for  higher  frequency 
signals,  and  correspondingly  lower  at  lower  frequencies,  with  the  maximum  ratio  being  8,  and  the  minimiun  being 
1.09.  This  example  was  at  appoximately  350  Mhz. 

We  would  like  to  study  the  behavior  of  HRR  and  our  precursor  based  hybrid  HRR  through  an  absorbing  and 
dispersive  media.  We  outline  the  assumptions  and  rules  that  will  be  used  throughout  this  comparison  below: 

•  The  noise  figure  stated  is  the  signal  to  noise  ratio  of  a  freespace  system  without  absorption  and  dispersion. 
This  recognizes  that  there  is  inherent  ambient  noise  in  the  radar  system,  and  that  the  medium  will  then  create 
a  much  lower  signal  to  noise  ratio. 
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Real  and  Imaginary  parts  of  the  Transmitted  Signals 


Real  and  Imaginary  parts  of  the  Received  Signals 


Figure  2:  We  present  the  signals  that  are  transmitted  for  the  precursor  HRR  above.  The  real  and  imaginary 
components  of  the  transmit  signal  are  presented  at  the  top.  We  swap  the  phase  of  the  signals  at  every  other  zero 
crossing.  Note  also  that  the  real  and  imaginary  components  of  the  transmit  signal  are  tt  out  of  phase  with  each 
other,  as  compared  to  7r/2  for  a  standard  exponential.  This  is  necessary  so  that  the  received  signals,  shown  at  the 
bottom,  will  be  7r/2  out  of  phase,  enabling  standard  HRR  processing  of  the  received  signals.  Note  that  the  precursor 
behavior  has  resulted  in  a  frequency  halving  of  the  signal,  but  the  signals  will  be  nearly  constant  phase  after  mixing 
with  a  propier  exponential  (  they  are  nearly  exponentials  themselves  ). 

Returned  HRR  signal 


Figure  3:  We  illustrate  the  effect  of  our  material  model  above.  A  standard  HRR  signal  was  transmitted,  at  approx- 
imatly  350  Mhz,  and  the  resulting  received  signal  is  displayed  at  the  top.  Note  the  strong  precursor  effect  at  the 
leading  and  trailing  edge.  We  believe  that  the  foliage  model  utilized  here  is  consistent  with  measmed  precursor  effects 
through  realistic  foliage.  Note  that  the  basis  of  the  precursor  HRR  is  to  induce  the  increased  power  shown  above  at 
the  leading  and  trailing  edges,  throughout  the  pulses.  This  is  shown  in  the  bottom  figure,  where  the  corresponding 
received  signal  firom  the  precursor  HRR  system  is  displayed.  There  is  5.8  times  more  power  in  the  precursor  signal, 
pving  hojie  for  increased  system  performance. 
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•  We  will  study  standard  HRR  and  our  precursor  HRR  with  various  base-band  ( lowest ),  frequency  settings.  We 
have  chosen  2D0,  300,  and  600  Mhz  as  our  initial  basebands,  with  approximately  500  Mhz  of  frequency  utilized 
above  each  of  these  basebands.  We  believe  that  this  is  a  realistic  comprimise  range  between  antennae  size  and 
power  considerations. 

•  In  each  study,  there  will  be  128  pulses  which  are  128  feet  long,  appropriately  spaced  in  frequency  via  HRR 
rules  to  yield  1  foot  of  resolution  in  a  freespace  HRR  system. 

•  There  are  two  gross-range  bins  included  in  our  study  (  128  feet  long  ),  which  are  then  resolved  into  128  fine 
range  bins  (  at  least  in  the  freespace  example  ).  Since  the  ground  pr^ents  a  defimte  base  for  range  gating,  we 
feel  this  is  more  than  sufficient. 

•  We  utilize  one  point  target,  placed  approximatly  31  feet  into  the  first  gross  range  bin. 

•  We  will  report  the  emperical  resolution  of  each  trial.  This  is  defined  for  our  purposes  as  being  the  number  of 
fine  range  bins  which  exceed  1/3  of  the  peak  value  of  the  transfer  function.  Thus  if  the  peak  value  is  150,  and 
there  are  11  fine  range  bin  values  around  the  i>eak  of  the  transfer  function  over  50,  the  resolution  would  be 
reported  as  being  11  feet. 

We  begin  this  comparison  by  studying  the  lowest  frequency  problem.  The  baseband  frequency  is  200  Mhz,  and 
approximately  2(K)-700  Mhz  of  frequency  range  is  utilized.  In  Figure  4,  we  consider  the  precursor-HRR  and  the 
standard  HRR.  Note  that  detection  for  this  point  target  is  possible  with  either  system,  but  that  the  precursor  HRR 
had  a  much  higher  emperical  resolution  of  3  feet,  as  opposed  to  11  feet  for  the  standard  HRR. 

To  try  and  close  this  resolution  gap,  we  utilized  the  standard  assumption  of  HRR,  that  equal  power  will  return 
from  each  frequency,  to  invert  the  spreading  of  the  transfer  function  and  try  to  reachieve  high  range  resolution.  The 
result  is  shown  in  Figure  5.  As  one  might  expect,  noise  from  the  higher  frequency  systems  has  essentially  eliminated 
the  abihty  to  locate  the  target.  Thus  we  were  not  able  to  get  the  standard  HRR  system  to  perform  at  the  same  noise 
levels,  with  the  same  resolution  as  the  precursor  system. 

Figure  6,  is  the  analogue  of  Figure  4,  with  the  baseband  being  300  Mhz.  Once  again  the  targets  are  identifiable, 
but  the  precursor  HRR  has  much  higher  emp>erical  resolution  (  3  feet  ),  than  the  standard  HRR  (11  feet  ).  Figure 
7  is  the  analogue  of  Figure  5,  at  300  Mhz.  We  attempted  to  recover  the  resolution  from  the  standard  HRR  but  the 
noise  levels  were  prohibitive. 

Figure  8  shows  a  comparisons  of  the  systems  with  a  base  frequency  of  600  Mhz,  and  is  analogous  to  Figures 
4  and  6.  Note,  however,  that  noise  has  absolutely  overcome  the  HRR  system,  making  target  identification  nearly 
impossible.  In  addition,  we  removed  the  noise  to  test  the  resolution  of  the  two  s3rstenis.  The  precursor  HRR  had  a 
3  foot  resolution,  just  as  at  200  and  300  MHz,  while  the  resolution  of  the  standard  HRR  system  had  dropped  to  22 
feet,  making  target  identification  nearly  impossible. 

We  illustrate  the  singular  values  of  the  systems  at  600  Mhz  baseband  in  Figure  9.  The  singular  values  of  the 
HRR  system  are  much  smaller  than  those  of  the  precursor  HRR  system,  explaining  the  large  differences  which  were 
noted  in  Figure  8.  Recall  that  larger  singular  values  imply  transmission  of  more  power  and  information. 
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Figure  4:  Above  we  present  a  comi>arison  of  the  precursor  HRR,  on  top,  vs.  standard  HRR  with  at  200  Mhz  below. 
While  both  systems  identify  the  targets,  the  precursor  HRR  has  an  emperical  resolution  of  3  feet,  while  the  standard 
HRR  has  an  emperical  resolution  of  11  feet. 


Noise  level  is  5  dB 


Figure  5:  Comparison  of  the  precursor  HRR,  vs.  standard  HRR  with  inversion  attempted  on  the  HRR,  at  200  Mhz. 
Note  that  the  standard  HRR,  with  frequency  equalization,  is  too  noisy  for  realistic  use. 
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Noise  level  is  5  dB 


Figure  6:  Above  we  present  a  coinp>arison  of  the  precursor  HRR  at  top,  vs.  standard  HRR  at  300  Mhz.  As  in  Figure 
4,  both  identify  the  target,  but  the  precursor  HRR  has  a  resolution  of  3  feet,  versus  a  resolution  of  11  feet  for  the 
standard  HRR. 


Noise  level  is  5  dB 


Figure  7:  Above  we  present  a  comparison  of  the  precursor  HRR  at  top,  vs.  standard  HRR  at  300  Mhz.  As  in  Figure 
5  we  attempted  to  recover  the  resolution  from  the  HRR  system,  and  noise  effects  made  this  non- viable. 


7 


Precursor-HRR,  Noise  level  5  dB 


Figure  8:  Above  we  present  a  comparison  of  the  precursor  HRR  at  top,  vs.  standard  HRR  with  at  600  Mhz.  As 
opposed  Figures  4  and  6,  the  standard  HRR  cannot  identify  the  target.  The  precursor  HRR,  however,  demonstrates 
robustness  and  still  has  an  emp)erical  resolution  of  3  feet.  Testing  the  resolution  of  the  HRR  system,  with  essentially 
no  noise,  revealed  that  it  was  22  feet.  Since  the  noise  is  already  prohibitive  in  the  HRR  system,  we  did  not  attempt 
to  recover  the  resolution  of  the  system  as  in  Figures  5  and  7. 
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Figure  9:  Comparison  of  the  singular  values  of  precursor  HRR,  in  blue,  and  traditional  HRR,  in  green,  with  the 
lowest  frequency  being  600  Mhz.  The  norm  of  the  singular  values  of  the  precursor  HRR  are  25  times,  i.e.  the 
Hilbert-Schmidt  norm  of  the  system,  is  25  times  larger  than  that  of  the  standard  HRR  system. 
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Noise  level  is  0  dB 


Figure  10:  Above  we  present  a  comparison  of  the  precursor  HRR,  on  top,  vs.  standard  HRR  with  at  200  Mhz  below. 
This  is  analogous  to  Figure  4,  but  with  a  noise  level  of  0  dB  instead  of  5  dB  in  Figure  4.  Note  that  the  HRR  system 
did  not  identify  the  target  at  all,  while  hte  precursor  HRR  system  accurately  performed  with  this  high  noise  level. 
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Noise  level  is  -5  dB 


Figure  11:  Above  we  present  the  output  of  the  precursor  HRR  at  200  Mhz,  with  a  very  low  signal  to  noise  ratio  of 
-5  dB.  Note  that  the  standard  HRR  system  could  not  accurately  obtain  the  resolution  obtained  here,  namely  3  feet, 
at  +5  dB.  The  precursor  HRR  system  ac<mrately  detects  the  system  at  -5  dB,  indicating  a  10  dB  win. 


Noise  level  is  0  dB 


Figure  12:  Above  we  present  the  output  of  the  precursor  HRR  at  600  Mhz,  with  a  low  signal  to  noise  ratio  of  0  dB. 
Note  that  the  standard  HRR  system  could  not  accurately  obtain  the  resolution  obtained  here,  namely,  at  +5  dB.  In 
addition  its  transfer  function  had  a  resolution  of  only  20  feet. 


10 


.  Derivation  of  Precursors  via  Finite  Toeplitz  Forms 

Tim  Olson 

Department  of  Mathematics 
The  University  of  Florida 
364  Little  Hall 
Gainesville,  FL  32611-8105 

ABSTRACT 

In  this  paper  we  derive  ’’precursors”  from  basic  principles  of  linear  algebra  and  operator  theory.  Evidence  is  given 
that  these  are  indeed  the  Brillouin  precursors.^  We  give  numerous  numerical  examples  and  illustrate  that  by  imposing 
a  finite  t.imp  limit  the  signals,  the  precursor  phenomena  will  be  generated  from  the  infinite  Toeplitz  form. 

1.  Introduction 

Precursors  were  first  explained  in  1914  by  Sommerfeld  and  Brillouin.  Their  primary  focus  was  on  analyzing 
these  pulses  with  resi)ect  to  causality.  The  conclusion  of  these  papers  was  that  the  leading  edge  of  the  Sommerfeld 
precursor  travels  at  the  vacuum  sf)eed  of  light,  but  is  absolutely  causal.  Similarly,  the  Brillouin  precursor  tends  to 
have  a  majority  of  its  energy  at  the  leading,  and  corresponding  trailing  edges  of  the  pulse  and  is  causal.  Unfortimately 
Brillouin  concluded,  due  to  incorrect  asymptotics,  that  these  pulses  would  not  be  useful.  Recent  work  has  shown 
that  Brillouin’s  analysis  was  flawed  and  that  the  Brillouin  precursor  is  the  dominant  term  in  the  far  field  of  any 
disp>ersive  and  absorbing  media.* 

The  methods  of  Brillouin,  Sommerfeld  and  Oughstun  were  asymptotic  in  nature.  In  that  work,  the  dominant 
nature  of  precursors  was  found  using  advanced  techmques  to  approximately  solve  a  complex  integral.  Often  in 
applied  mathematics,  we  strive  to  understand  the  structure  of  transmission  operators  via  operator  theory.  If  the 
transmission  operator  can  be  examined  as  a  compact  operator,  the  generated  structure  is  very  informative.  The 
inputs  and  outputs  to  the  operator  are  separated  into  orthogonal  subspaces,  with  the  power  passing  through  each 
subspace  clearly  described  by  the  singular  values,  which  are  analagous  to  eigenvalues.  We  will  show  that  Brillouin 
precursors,  or  pulses  remarkably  similar  to  Brillouin  precursors,  are  the  dominant  singular  vectors  associated  with 
transmission  through  media  such  as  foliage  or  water  using  basic  techniques  of  operator  theory  or  linear  algebra. 

2.  Background  material  and  notation 

We  begin  in  the  spirit  of  Slepian  and  Pollack.*  In  that  work,  the  singular  functions  of  the  joint  time-frequency  cut 
off  operators  were  derived,  resulting  in  the  prolate  spheroidal  wave  functions.  We  will  consider  operators  which  are 
similar,  but  more  general  than  the  joint  time-frequency  operator.  These  operators  are  physically  motivated  by  a  wide 
variety  of  electromagnetic  propagation  problems.  The  convolution  operators  which  describe  the  evolution  of  a  pulse 
r(t,  z)  through  a  homogeneous  linear  medium  have  a  very  simple  form.  Given  an  initial  plane  wave  signal  which  is 
incident  on  a  homogeneous  mediiun,  s(<,0)  =  s(t),  the  pulse  at  time  t,  and  distance  z  is  appropriately  modeled  by 

r(t,  ^)  =  j  s(t)Az(t  -  t)dT  =  Lz(s(t)).  (1) 

Unless  otherwise  noted,  all  integrals  are  over  the  real  line.  Convolution  operators  L  of  the  type  (1)  have  been  heavily 
studied  and  are  well  understood.  The  Fourier  transform  diagonalizes  the  operator  and  the  spectrum  of  the  operator 
is  the  continuous  Fomier  transform  of  Az,  for  any  fixed  distance  z.  A  monochromatic  signal  s(t)  transmitted  at  a 
frequency  wq,  will  be  absorbed  according  the  real  part  of  the  Fourier  transform  Az(wo).  Dispersion  is  described  by 
the  complex  portion  of  .4z(w;o)-  When  the  signal  is  not  monochromatic,  then  the  resulting  signal  r{t,  z)  =  rz{t)  has 
a  Fourier  transform  which  is  the  product  of  A,  and  s,  i.e.  rz{w)  =  Az{w)s{w). 


This  work  was  supported  by  the  Air  Force  Office  of  Scientific  Research, 


Slepian  and  Pollack  utilized  the  finite  length  of  signals  to  alter  an  operator  of  the  type  (1).  This  alteration 
creates  a  compact  operator,  with  a  corresponding  discrete  set  of  singular  values  and  singular  vectors  as  opposed  to 
the  continuous  spectrum  of  L*.  Followingthis  development  we  will  consider  pulses  of  finite  length  I,  and  corresponding 
new  operators  L*,,  which  describe  a  pulse  evolving  through  a  distance  z  of  a  medium.  Formally  we  have 

/OO  /‘OO 

s{r)Az{T  -  t)dr  =  I  s{r)Xi{T)Az{r  -  t)dr  L],{s(t)),  (2) 

-OO  Joo 

where  Xi{t)  =  1  for  r  e  [0,  /],  and  is  0  otherwise.  Thus  our  old  kernel  was  Az{t)  and  our  new  kernel  is  t)  = 
Xi  {t)Az  (T-t).  A  basic  result  of  functional  analysis  states  that  when  a  kernel  of  the  type  If  ‘  (t ,  t)  is  square  integrable, 
the  corresponding  operator  is  a  compact  operator.  This  is  stated  clearly  in 

Theorem  1  (The  HilbERT-SchMIDT  Theorem;).  Let  an  operator  L  be  defined  by 

Li/m  =  I  f{T)Git,r)dr  (3) 

and  let  \\G{t,T)\\2  <  oo.  Then  L  is  a  compact  operator,  and  it  follows  that  there  exist  singular  vectors  and  singular 
values 

/oo 

/(r)G(t,  T)dT  =  '^<rk  if,  Vk)uk.  (4) 


The  values  au  are  called  the  singular  values  and  the  vectors  {wfe}  and  {ufc}  are  correspodingly  called  the  left  and  right 
singular  vectors.  Moreover,  the  energy  of  the  singular  values  is  exactly  that  of  the  kernel,  or 


[  f  \Git,T)fdtdr  = 


(5) 


At  this  time,  let  us  adopt  some  notation.  We  are  dealing  with  a  class  of  compact  operators  which  is  indexed 
by  the  initial  finite  signal  length  I,  and  the  propagation  distance  z,  so  we  refer  to  the  kernels  of  the^  o^rators  as 
if' .  Similarly  we  will  refer  to  the  corresponding  singular  values  (Tfc(z,  1),  where  k  is  the  index,  and  likewise  the  left 
angular  vectors  «fc(z,  1),  and  right  singular  vectors  Vk{z,  1).  Thus  k  runs  from  0  to  oo,  and  the  necessarily  positive 
angular  values  decrease  by  convention.  Thus  the  dominant  left  and  right  singular  vectors  are  always  uoik,l),  and 
«o(fe,  ly  We  will  also  refer  to  the  transmission  operator,  without  regard  to  the  finite  pulse  length  as  L*. 

.3.  Basic  Results: 

With  basic  development  of  theory,  we  will  now  calculate  a  specific  example  from  third  party  electronaagnetic  data. 
The  data  we  have  utilized  was  provided  by  the  computational  electromagnetics  group  at  Brooks  Air  Force  b^, 
and  this  specific  data  set  is  attributed  to  Grant.  The  transfer  function  generated  by  the  absorption  coefficient,  i.e. 
exp{—a(w))  is  shown  in  Figme  3. 

We  will  now  analyze  the  Angular  value  decompositions  of  the  generated  operators  at  various  distances,  and  with 
various  initial  pulse  widths.  To  begin  with  we  start  out  with  a  pulse  which  is  only  7.5  cm  long.  Thus  we  will  be 
considering  the  optimal  pulse  to  transmit  one  meter  through  water,  assuming  that  the  initial  pulse  is  less  than  7.5 
cm  long.  The  spectrum  and  time  plot  of  this  dominant  singular  vector,  uo{7.5, 1),  is  shown  in  Figure  3. 

We  have  arbitrarily  picked  1  meter  as  our  propagation  distance,  and  7.5  centemeters  as  our  initial  allowed  pulse 
length.  Natural  questions  at  this  time  revolve  around  these  arbitrary  choices,  namely: 

•  What  effect  does  the  initial  finite  pulse  length  have  on  the  singular  vectors  in  the  far  field? 

•  Does  the  dominant  singular  vector  at  1  meter  propagate  well  over  20  meters? 

We  will  suggest  answers  to  these  questions  now.  The  answers  seems  to  be  as  following 

•  Over  short  distances,  the  dominant  pulse  depends  on  length  and  propagation  distances. 

•  Over  longer  distances,  the  dominant  pulse  becomes  indej)endent  of  initial  pidse  length  and  propagation  distance. 

Let  us  look  at  a  couple  of  examples  to  support  these  claims. 
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Figure  1:  Testing  the  water 
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Figure  2:  On  the  left  we  have  displayed  a  time  plot  of  the  dominant  angular  function,  and  the  second  singular 
function.  Note  that  the  dominant  sin^ar  (  in  blue  )  function  looks  very  similar  to  the  leading  edge  of  a  Brillouin 
precursor.  On  the  right  we  have  displayed  the  transfer  function  for  the  operator  ( in  blue  ),  and  the  absolute  value 
of  the  spectrum  of  the  dominant  singular  function  and  the  second  singular  function  (  in  green  and  red  respectively 
).  Note  that  the  spectrum  of  the  dominant  singular  function  is  very  nearly  that  of  the  transfer  function  or  kernel. 

3.1.  How  does  propagation  distance  affect  the  singular  vectors? 

To  begin  with,  we  want  to  compare  the  propagation  of  two  different  singular  vectors;  1)  The  dominant  singular  vector 
for  propagating  through  1  meter  of  water,  and  2)  the  dominant  singular  vector  for  propagating  through  20  meters 
of  water.  We  will  fix  the  pulse  length  at  7.5  centimeters. 

The  test  is  set  up  in  the  following  way.  We  calculate  the  singular  vectors  in  the  above  to  situations.  Then  we 
propagate  them  through  a)  1  meter  of  water,  and  b)  50  meters  of  water.  The  first  singular  vector  was  calculated 
to  be  optimal  for  experiment  a)  or  1  meter  of  water,  and  the  second  singular  vector  for  propagation  through  b)  50 
meters  of  water.  The  question  is:  ”Do  they  both  perform  well?”. 

The  answer  was  very  conclusive.  The  results  were  nearly  identical  after  propagating  through  the  various  dis- 
tanceds.  This  can  be  seen  in  Figures  3  and  4  where  we  tested  the  results  through  1  and  20  meters  of  water 
respectively.  Note  that  the  time  plots  of  these  singular  vectors  are  nearly  indistinguishable.  The  frequency  plots 
show  the  transfer  function,  the  first  singular  vector,  and  the  second  singular  vector  for  each  case.  Note  that  the 
dominant  ”  precursor”  pulse  in  eadi  case  is  indistinguishable  whether  it  was  originally  calculated  at  1  meter,  or  50 
meters,  but  the  second  singular  vectors  are  orthogonal  and  therefore  very  different. 

We  then  asked  the  question:  ”How  close  are  the  singular  vectors  to  being  identical  when  they  are  calculated  at 
different  distances,  originally,  and  then  propagated  toward  the  far  field?”  The  answer  to  this  is  given  in  Figure  5.  The 
singular  vectors  calculated  at  a  short  distance,  and  a  very  long  distance,  are  extremely  similar,  with  all  correlations 
>  .99. 

The  conclurion  of  this  experiment  is  that  the  dominant  singular  vectors,  calculated  for  propagation  distances, 
projjagate  to  become  nearly  the  dominant  singular  vectors  at  the  next  distance.  Thus  they  appear  to  be  precursors. 
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Figure  3;  On  the  left  we  have  displayed  a  time  plot  of  the  dominant  singular  function  after  traveling  1  m  through 
water  (  in  blue  ),  and  the  dominant  singular  function  for  propagating  through  20  m,  after  it  has  propagated  through 
1  m  of  water.  Symbolically  this  is  «o(7.5, 1)  =  L^(uo(7.5, 1)),  compared  against  a  second  vector  L^{yo{7.^,  20)).  Note 
that  the  dominant  singular  calculate  for  20  m  of  water,  is  also  nearly  equal  to  the  dominant  singular  vector  for 
propagation  through  1  m  of  water.  Thus  the  arbitrary  dioice  of  propagation  distance  does  not  seem  to  matter.  The 
dominant  singular  vector  for  propagating  through  20  m,  had  a  correlation  with  the  other  singular  vector  of  .997. 
The  red  curve  on  the  right  is  that  of  the  second  singular  function.  Since  it  is  orthogonal  to  the  first  singular  function 
calculated  at  1  m,  it  will  be  nearly  orthogonal  to  the  singular  function  calculated  at  20  m,  because  of  the  noted 
correlation  above. 


Propasotlon  <liManc«  20  m  Fraqyency  Plot 


Figure  4:  On  the  left  we  have  displayed  a  time  plot  of  the  dominant  singular  function  after  traveling  20  m  through 
water  (  in  blue  ),  and  the  dominant  singular  function  for  propagating  through  1  m,  after  it  has  propagated  through 
20  m  of  water.  Symbolically  this  is  wo(7.5, 20)  =  L“(t;o(7.5, 20)),  compared  against  a  second  vector  I,^°(wo(7.5, 1)). 
Note  that  the  dominant  singular  vector  calculated  for  1  m  of  water,  is  also  nearly  equal  to  the  dominant  singular 
vector  for  propagation  through  20  m  of  water.  Thus  the  arbitrary  choice  of  propagation  distance  does  not  seem 
to  matter.  The  dominant  singular  vector  for  propagating  through  20  m,  had  a  correlation  with  the  other  singular 
vector  of  .997.  Once  again  the  red  curve  at  the  right  is  the  second  singular  function  calculated  at  20  m,  which  will 
be  orthogonal  to  the  dominant  singular  functions. 
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Correlations 


Figure  5:  Above  we  Ulustrate  the  projection  of  the  dominant  left  singular  vector  calculated  at  1  m,  onto  the 
dominant  left  singular  vectors  for  distances  1-50.  Symbolically  we  are  calculating  the  projection  of  L*’(vo(7.5, 1)) 
onto  uo(7.5,  k)  =  L*'(t;o(7.5,  k)).  Note  that  there  is  almost  no  difference  in  the  angular  vectors. 


Figure  6:  Above  we  illustrate  the  projection  of  the  dominant  left  singular  vector  calculated  at  distances  1-50  m, 
with  initial  time  contraints  of  either  7.5  cm,  or  15  cm.  Symbolically  we  are  calculating  the  projection  of  mo(7.5,  k)  = 
L*'(tto(7.5,  k))  onto  mo(15,  k)  =  L*'(uo(15,  k)).  The  initial  pulse  length  seems  to  become  irrelevant  in  the  far  field,  as 
the  distance  increases,  as  is  illustrated  by  the  correlation  gomg  to  one,  i.e.  the  left  smgular  vectors  are  converging 
to  one  another. 

Our  study  of  the  dependence  of  the  dominant  singular  vectors  on  distance,  has  indicated  that  the  dominant 
subspace  which  is  generated  by  the  first  singular  vector  seems  to  propagate.  In  other  words  if  we  laimch  a  right 
gingnlar  vector  into  the  material,  the  output  at  1  meter  is  dominant  left  singular  vector.  Very  importantly,  however, 
the  results  of  the  last  section  suggest  that  if  we  let  this  vector  propagate  an  additional  19  meters,  it  will  still  be  nearly 
the  dominant  singular  vector  for  the  medium  at  20  meters.  Thus  we  do  not  have  to  know  the  nature  of  the  medium, 
or  perhaps  the  depth  of  the  foliage,  exactly  in  order  to  lauch  a  pulse  which  is  close  to  optimal  from  a  penetration 
viewpoint. 

3.2.  Dependence  on  the  initial  time  length 

We  now  study  the  dependence  of  the  singular  vectors  on  the  initial  pulse  length.  Since  we  arbitrarily  chose  the  initial 
pulse  length  of  7.5  cm,  we  need  to  imderstand  how  this  initial  choice  affects  the  calculated  singular  vectors. 

To  do  this,  we  calculated  the  optimal  singular  vectors  first  at  a  distance  z,  with  initial  pulse  length  7.5  cm, 
and  then  at  distance  z,  with  initial  pulse  length  15  cm.  These  results  are  shown  in  Figure  5.  Note  that  since  the 
corr^ation  converges  towards  1,  it  is  safe  to  assume  that  the  singular  vectors  are  converging  to  a  single  precursor  in 
the  far  field. 
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Absorption  coefficient  alpha(>iy) 


Figure  7:  Above  we  illustrate  the  absorption  coefficient  which  was  utilized  for  the  computations,  as  a  function  of 
frequency.  Notice  that  it  is  nearly  linearly  dependent  upon  frequency 


4.  Asymptotic  Results 

We  will  now  prove  one  of  the  major  theorems  of  this  paper.  Much  of  the  fascination  with  precursors  is  due  to  the 
fact  that  they  propagate  with  an  absorption  which  is  rather  than  We  will  now  give  a  simple  proof  of 

this  fact,  under  a  basic  h3q>othesis  which  seems  valid  in  the  case  of  the  water  data  which  we  have  just  presented. 

Our  primary  assumption,  for  the  purpose  of  this  r^ult  is  that  the  absorption  coefficient  Qf(w7)  is  linearly  dep)endent 
upon  the  frequency  for  low  frequencies,  or  ^  kw  for  w  small  where  fe  is  a  positive  constant.  We  now  present 
the  following 

Theorem  2  (Asymptotic  Decay).  Assume  that  propagation  within  a  dielectric  material  is  modeled  correctly 
by  (2).  Assume  further  that  the  absorption  coefficient  a{w)  kw,  where  k  is  a  positive  constant  Then  the  energy 
of  the  kernels  associated  vnth  propagation  of  a  distance  z,  K[{t,  r),  decays  at  the  rate  0(2:“^/^).  Moreover,  by  the 
HUbert-Schnidt  Theorem,  it  follows  therefore  that  the  sum  of  the  squares  of  the  singular  values  also  decays  at  a  rate 
ofO(z-^l’^). 

Proof.  The  proof  is  very  straight  forward  manipulation  of  the  integrals,  since 


^00  /»oo  poo  POO 

/  /  \Kl(t,T)fdtdT  =  /  /  \Xi(t)Az{T-t)\^dtdT 

J  ~oo  J—00  J—00  00 

=  /  /  \Xi{t)A:,{T -t)\^dtdr 

J-00  J-i 

=  f  f  \Az{t —  t)\^dTdt 

J-l  J-00 

-  f  f  \Az{T)\‘^dTdt 

J-iJ-00 

/OO 

\Az{r)\^dT. 

-OO 


(6) 


Now  the  Fourier  Isometry  and  (6)  gives  us  that 
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Figure  8:  We  illustrate  the  asymptotics  of  the  singular  vector  which  we  believe  to  be  the  precursor  field  above.  We 
have  plotted  ^{z)  *  .25  +  1.5  in  blue,  and  l/a(7.5,  z)  in  green.  Note  that  the  fit  is  very  good,  emphasizing  that  the 
gingiilnr  vector,  which  we  believe  to  be  the  precursor,  does  indeed  display  0{z~'^l^)  behavior. 

=  4/  r  (7) 

Jo 

From  the  substitution  u  =  kwz,  it  follows  that  dw  =  1/ (kz)du,  and 

\\Kif  =  4l^j^e-^^du  =  0{z-^).  (8) 

By  taking  a  square  root,  we  end  up  with  the  result  that  |lii:'||  =  which  corresponds  with  the  result  of 

Brillouin.  □ 

While  we  have  shown  that  the  rate  of  decay  of  the  operator  is  0(z~^^^),  we  have  not  yet  proven  that  the  dominant 
singular  values  decay  at  this  rate.  It  is  the  case  that  each  frequency  in  the  kernel  decays  exponentially,  except  for 
the  origin,  but  the  energy  of  the  kernel  decays  much  more  slowly,  i.e.  0{z~^^^).  Prom  the  Hilbert-Schmidt  Theorem 
we  know  that  the  sum  of  the  squares  of  the  singular  values  decays  as  0(2:“^/^),  but  perhaps  each  of  them  decay 
exponentially  just  as  the  frequencies  in  the  oi>erator. 

To  study  this,  we  present  a  final  numerical  computation.  We  begin  by  showing  that  the  singular  vector  associated 
with  the  precursor  does  indeed  display  aymptotics  which  seem  to  be  0{z~^^^).  This  is  seen  in  Figure  8.  Secondly, 
we  will  compare  the  amount  of  energy  in  the  precursor  field,  to  the  total  energy  in  the  operator,  which  is  seen  in 
Figure  9. 

5.  Relationship  to  Standard  Brillouin  Precursors 

We  have  shown  that  basic  operator  theory  generates  the  dominant  singular  vector  for  transmission  through  a  medium, 
and  that  this  dominant  singular  vector  looks  qualitatively  similar  to  a  Brillouin  precursor.  We  have  also  shown  that 
this  dominant  singular  vector  seems  to  be  relatively  robust  to  changes  in  the  finite  time  restriction,  and  distance 
through  the  material. 

We  have  not  directly  compared  this  dominant  angular  vector,  which  we  believe  to  be  the  Brillouin  precursor,  to  a 
Brillouin  precursor  computed  in  a  traditional  manner.  We  present  this  experiment  below.  In  Figure  10  we  illustrate 
the  result  of  propagating  a  square  windowed  sinusoid  through  an  absorbing  and  dispersive  media.  Note  the  two 
Brillouin  precursors  at  the  leading  and  trailing  edges.  In  Figure  11  we  compare  the  Brillouin  precursor,  generated  as 
in  Figure  10,  to  the  dominant  singular  vector  for  the  transmission  operator,  which  we  claim  is  the  Brillouin  precursor 
or  nearly  the  Brillouin  precursor.  Note  that  they  are  very  similar. 

We  will  now  examine  whether  the  dominant  singular  vector  becomes  more  or  less  similar  to  the  Brillouin  precur¬ 
sor,  as  distance  increases.  The  question  being  ”Are  they  asymptotically  the  same”.  The  literature  provides  great 
motivation  for  this  question.  It  has  been  proven  that  the  Brillouin  precursor  is  the  dominant  effect  in  the  far  field. 
The  dominant  singular  vector  is  by  definition  the  dominant  effect.  Science  leaves  little  room  for  these  two  dominant 
effects  to  be  dissimilar.  The  experiment  illustrated  in  Figure  12  suggests  that  yes,  they  are  the  same  in  the  far  field. 
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Energy  Distribution  in  the  Rrocursor  Field 


Figure  9:  We  have  tested  the  dominance  of  the  precursor  field  in  the  above  diagram.  We  have  plotted  the  ratio, 
111^112/^2(7.5, 0),  or  the  ratio  of  the  operator  norm  to  the  value  of  the  dominant  singular  value.  Note  that  this  ratio 
is  strictly  bounded  above  by  1,  and  that  the  ratio  rapidly  approaches  1,  starting  out  at  over  .9.  Thus  nearly  all  of 
the  energy  of  the  operator  is  quickly  in  the  precursor  field. 


RscatvBd  Signal 


Figure  10:  In  this  figure  we  illustrate  two  classical  Brillouin  precursors,  the  first  at  the  beginning  of  the  pulse,  and  the 
second,  with  a  negative  orientation  at  the  end  of  the  pulse.  This  is  typical  of  the  effect  realized  when  one  transmits 
a  square  windowed  sinusoid  through  and  absorbing  and  dispersive  media. 


Precursor  Comparisons 


Figure  11;  In  this  figure  we  illustrate  the  singular  vectors  of  the  transmission  operator,  which  we  believe  to  be 
essentially  Brillouin  precursors,  against  a  Brilluin  precursor  which  was  generated  by  sending  a  square  wave  through 
the  same  finite  distance  of  an  absorbing  media,  such  as  in  Figure  10.  Note  that  they  are  extremely  close  in  size  and 
shape. 
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10  meters  20  meters 


Figure  12:  Above  left  we  illustrate  the  singular  vectors  of  the  transmission  operator  at  10  meters,  and  Brillouin 
precm^r  which  was  generated  by  sending  a  square  wave  through  10  meters  of  an  absorbing  media.  Above  right 
we  illustrate  the  same  thing,  but  after  20  meters  of  material  Note  that  they  are  originally  close,  but  are  becoming 
closer  with  distance  suggesting  asymptotic  convergence. 

6.  Conclusion: 
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Abstract 

A  novel  time  domain  solver  of  Maxwell’s  equations  in  passive  (dispersive  and  ab¬ 
sorbing)  media  is  proposed.  The  method  is  based  on  the  path  integral  formalism  of 
quantum  theory  and  entails  the  use  of  (z)  the  Hamiltonian  formalism  and  (iz)  pseu- 
dospectral  methods  (the  fast  Fourier  transform,  in  particular)  of  solving  differential 
equations.  In  contrast  to  finite  differencing  schemes,  the  path  integral  based  algorithm 
has  no  artificial  numerical  dispersion  (dispersive  errors),  operates  at  the  Nyquist  limit 
(two  grid  points  per  shortest  wavelength  in  the  wavepacket)  and  exhibits  an  expo¬ 
nential  convergence  as  the  grid  size  increases,  which,  in  turn,  should  lead  to  a  higher 
accuracy.  The  Gauss  law  holds  exactly  with  no  extra  computational  cost.  Each  time 
step  requires  0{N  log2  ^)  elementary  operations  where  N  is  the  grid  size.  It  can  also 
be  applied  to  simulations  of  electromagnetic  waves  in  passive  media  whose  proper¬ 
ties  are  time  dependent  when  conventional  stationary  (scattering  matrix)  methods  are 
inapplicable.  The  stability  and  accuracy  of  the  algorithm  are  investigated  in  detail. 


^electronic  mail:  shabanovOphys. ufl.edu 


1 


1 


Introduction 


In  this  study  a  time  domain  solver  of  Maxwell’s  equation  in  passive  (dispersive  and  absorbing) 
media  is  developed.  The  main  motivation  of  this  work  is  to  bring  methods  of  computational 
quantum  physics  into  classical  electromagnetic  theory.  One  of  the  great  advantages  of  time 
domain  methods  over  stationary  (scattering  matrix)  methods  is  that  a  single  simulation  of 
the  scattering  of  a  wide  band  wave  packet  can  determine  basic  physical  properties  of  the 
target  (e.g.,  transmission  and  reflection  coeflicients)  in  the  entire  frequency  band  covered  by 
the  initial  wave  packet.  Time  domain  methods  also  allow  for  a  unique  possibility  to  observe 
all  immediate  effects  on  fields  caused  by  the  target  or  by  a  surrounding  passive  medium, 
which  greatly  facilitates  qualitative  understanding  of  the  interaction  of  an  electromagnetic 
pulse  with  media  and  targets.  Another  important  advantage  is  that  the  target  geometry 
(or  medium  ph3^ical  properties)  may  vary  with  time  and  this  time  dependence  cannot  be 
removed  by  going  over  to  a  moving  reference  frame.  Stationary  methods  are  simply  inappli¬ 
cable  to  these  kind  of  problems. 

Erom  the  computational  point  of  view,  the  proposed  approach  is  based  on  pseudospectral 
methods.  The  essential  advantages  of  pseudospectral  algorithms  over  conventional  finite  ele¬ 
ment  or  finite  difference  schemes  in  solving  differentird  equations  are  [1]:  (i)  the  exponential 
versus  polynomial  rate  of  convergence  as  the  grid  size  (or  the  basis  dimension)  increases;  (it) 
the  absence  of  dispersive  errors  and  (iti)  efficiency  in  numerical  calculations.  Time  domain 
algorithms  in  combination  with  pseudospectral  methods  have  become  the  state-of-the-art 
technique  in  numerical  studies  of  quantum  dynamics  by  solving  the  corresponding  initial 
value  problem  for  the  Schrodinger  equation  (see,  e.g.,  [2]).  A  typical  sdgorithm  entails  an 
approximate  computation  of  an  object  called  the  path  integral  (or  functional  integral)  intro¬ 
duced  by  Feynman  [4].  Here  Maxwell’s  theory  in  general  dispersive  media  is  reformulated 
in  the  Schrodinger  (Hamiltonian)  formalism.  Then  the  path  integral  formalism  is  applied 
to  the  initial  value  problem  in  Maxwell’s  theory  in  passive  media  to  develop  a  numerical 
algorithm.  The  main  objective  of  this  work  is  to  give  a  theoretical  assessment  of  the  path 
integral  based  solver  of  the  initial  value  problem  for  Maxwell’s  equations.  Numerical  tests 
and  applications  will  be  discussed  elsewhere  [3]. 

It  is  shown  here  that  basic  principles  of  the  path  integral  formalism  lead  to  a  true  time 
domain  algorithm  which  indeed  enjoys  the  advantages  of  pseudospectral  methods.  In  par¬ 
ticular,  among  the  aforementioned  features,  (*)  is  provided  by  the  use  of  the  fast  Fourier 
method  [5]  as  a  part  of  the  algorithm,  (it)  is  a  consequence  of  Nelson’s  construction  [6] 
of  the  path  integral  which  is  embedded  in  oxur  algorithm,  (tit)  is  due  to  the  fast  Foiuier 
method  and  some  analytical  results  that  speed  up  numerical  computations.  The  algorithm 
has  another  great  advantage  over  finite  difference  schemes:  The  Gauss  law  is  implemented 
exactly  with  no  extra  computational  cost  (Theorem  8.2).  For  widely  used  multi-resonant 
Lorentz  models  of  passive  media,  the  algorithm  is  unitary,  meaning  that,  the  energy  of  a 
wave  packet  is  preserved  exactly  in  dispersive  media  with  no  attenuation  (Theorem  6.1).  It  is 
also  unconditionally  stable  (Theorems  7.1  and  7.3)  versus  conditionally  stable  finite  element 
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or  finite  difference  algorithms  [1]  (see  also  [7]).  A  possible  drawback  of  the  algorithm  (to  be 
tested  numerically)  is  that  the  use  of  the  fast  Fourier  method  in  combination  with  Nelson’s 
construction  of  the  path  integral  might  require  additional  computational  costs  for  boundary 
value  problems  with  complicated  boundary  geometry.  In  our  approach,  conventional  boimd- 
ary  conditions  axe  not  imposed  on  electromagnetic  fields.  Targets  and  medium  interfaces  are 
modeled  by  discontinuous  medium  parameters.  The  problem  arises  firom  well  known  features 
of  the  Fourier  method  [5]:  Aliasing  and  low  convergence  rates  for  non-smooth  functions.  In 
this  study  we  oflfer  one  possible  way  to  deal  with  this  problem  while  keeping  the  Fourier 
basis  in  the  algorithm.  Alternative  pseudospectral  approaches  to  circumvent  the  problem 
exist  and  are  mentioned  here,  but  not  discussed  in  detail. 

The  basic  idea  of  the  path  integral  approach  to  solving  linear,  homogeneous,  evolutionary 
differential  equations  (numerically  or  analytically)  is  based  on  the  Hamiltonian  formalism. 
In  the  framework  of  the  Hamiltonian  formalism,  an  original  system  of  differential  equations 
is  transformed  to  an  equivalent  system  of  first-order  (in  time)  differential  equations  by  ex¬ 
panding  the  original  configuration  space,  that  is,  by  going  over  to  a  generalized  phase  space 
where  all  time  derivatives,  save  for  the  one  of  highest  order,  become  independent  variables 
[8].  A  generic  linear  homogeneous  first-order  system  can  be  written  in  the  form 

=  ’*'t=o  =  ^o,  (1-1) 

where  dt  stands  for  the  partial  derivative  with  respect  to  time  t,  a  linear  operator  Ti.  is  called 
Hamiltonian,  while  "3?^  is  called  a  state  vector  (or  wave  function).  It  is  an  element  of  the 
generalized  phase  space  of  the  system  and  viewed  as  a  collection  (colmnn)  of  the  original 
variables  and  their  time  derivatives.  The  generalized  phase  space  is  equipped  with  an  iimer 
product  and  becomes  a  Hilbert  space.  State  vectors  are  typically  vector-valued  functions  in 
R®,  and  the  Hamiltonian  is  a  differential  operator.  The  choice  of  the  inner  product  depends 
on  the  problem  at  hand.  One  usually  requires  that  componets  of  are  elements  of  ]L2(R®). 

In  general,  upon  going  over  to  the  Hamiltonian  formalism,  there  might  occur  constraints 
[9,  10] 

=  (1.2) 

with  Co  being  a  set  of  linear  operators  which  do  not  contain  time  derivatives;  a  enumerates 
the  constraint  operators.  The  constraints  must  be  preserved  in  the  time  evolution  which 
is  described  by  (1.1).  In  other  words,  the  solution  is  sought  in  the  subspace  of  the  Hilbert 
space  defined  by  (1.2).  Depending  on  the  type  of  constraints,  there  arc  different  ways  of 
developing  the  corresponding  path  integral  formafism.  In  Maxwell’s  theory,  the  constraint  is 
the  Gauss  law,  and  it  is  of  the  “first  class”  in  the  Dirac  terminology  [10].  The  characteristic 
feature  of  a  first  class  constrained  system  is  that 

[H,C„]~Ca,  [Ca,Ct]^Cc.  (1.3) 

A  consequence  of  (1.3)  is  that  if  the  initial  configuration  ’I'o  satisfies  the  constraints,  then  so 
does  the  solution  of  (1.1).  However,  after  the  projection  of  the  Hilbert  space  spanned  by 
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onto  a  finite-dimensional  subspace  (e.g.,  a  projection  on  a  subspace  associated  with  a  finite 
spatial  grid  as  is  done  in  Section  3),  which  is  required  for  numerical  simulations,  the  involu¬ 
tion  condition  (1-3)  can  be  violated  causing  problems  in  simulations.  For  instance,  the  Gauss 
law  is  typically  violated  in  any  finite  differencing  approach  to  simulations  of  electromagnetic 
wave  packet  propagation.  Special  efforts  have  to  be  made  to  ensure  the  transversality  of 
the  radiation  field  in  Maxwell’s  theory,  which,  in  turn,  comphcates  simulation  algorithms 
and  increases  computational  costs  (e.g.,  when  enforcing  the  Gauss  law  in  finite  difference 
schemes  on  the  grid  via  the  Lagrange  multiplier  method) .  It  is  one  of  the  advantages  of  the 
proposed  path  integral  based  algorithm  that  the  Gauss  law  can  be  strictly  enforced  with  no 
additional  computational  costs  for  generic  passive  media  (Theorem  8.2). 

The  solution  to  Eq.  (1.1)  is 

—  exp  {tH)  'Fo  =  Ut^o  ,  t  >  0  ,  (1-4) 

assuming  that  the  exponential  of  H  exists.  If  the  Hamiltonian  is  time  dependent  then  the 
following  replacement  has  to  be  made  in  (1.4) 

exp  {tH)  — >  T  exp  drH^  =  Ut  ,  (1-5) 

where  Texp  stands  for  the  time-ordered  exponential.  The  operator  Ut  is  defined  as  the 
fundamental  solution  of  (1.1),  dtUt  —  H^t  with  Ut=Q  being  the  identity  operator.  The  fun¬ 
damental  solution  has  the  semigroup  property,  Ut^+t^  =  Ut^t^>  The  action  of  the  evolution 
operator  Ut  on  the  initial  configuration  can  be  written  via  its  integral  kernel, 

'Ft(r)  =  /  dr'Wi(r,rO’Fo(r')  .  (1-6) 

Using  the  semigroup  property  of  the  evolution  operator,  the  entire  time  evolution  can  be 
viewed  as  consecutive  actions  of  the  infinitesimal  evolution  operator  Wai,  where  Af  is  a  time 
step.  If  the  kernel  of  the  infinitesimal  evolution  operator  is  known,  then  the  kernel  of  the 
evolution  operator  can  be  computed  as  the  convolution 

Wf(r, t')=  f  dri  •  •  •  dr„ ^At (r, r„)WAt(rn, **n-i)  •  •  'Um{^u ^')  (1-7) 

with  A*(n  +  1)  =  t.  The  integration  variables  can  be  regarded  as  points  r*  =  r(4),  where 
tk  =  kAt,  A:  =  0, 1, ...,  n-t- 1,  on  a  path  r(r)  connecting  points  r(T  =  f)  =  r  and  r(r  =  0)  =  r'. 
In  the  limit  At  0  the  convolution  (1.7)  can  be  viewed  as  a  sum  over  all  paths  coimecting 
the  initial  and  final  points.  This  is  the  gist  of  the  Feynman  path  integral  representation 
of  the  fundamental  solution  of  (1.1).  A  nontrivial  problem  is  to  find  the  measure  on  the 
space  of  paths.  For  example,  vfH  —  A  (the  Laplace  operator),  it  can  be  shown  that  the 
limit  exists,  and  that  the  measure  coincides  with  the  Wiener  measure  which  has  support  in 
the  space  of  all  continuous,  but  nowhere  differentiable  paths  (trajectories  of  the  Brownian 
motion)  pinned  at  the  end  points.  In  quantum  mechanics,  the  problem  is  more  subtle,  but 
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can  still  be  solved  [11].  The  existence  of  the  proper  measure  on  the  space  of  paths  opens 
up  an  attractive  possibility  to  use  Monte-Carlo  methods  of  computing  the  sum  over  paths 
which  is  the  gold  standard  algorithm  in  quantum  and  statistical  physics. 

However,  the  present  study  does  not  intend  to  tackle  the  measure  problem  for  the  path 
integral  representation  of  Maxwell’s  theory,  but  rather  offers  a  solution  of  a  more  modest 
problem.  Namely,  how  the  conventional  way,  outlined  above,  of  deriving  the  path  integral 
from  the  original  differential  equation  can  be  used  to  obtain  an  algorithm  for  numerical 
simulations  of  the  convolution  (1.7)  for  a  small,  but  finite  At.  Similar  ideas  have  been  applied 
to  non-dispersive  and/or  random  media  as  well  as  to  scattering  problems  and  waveguides  [12]. 
Our  approach  applies  to  general  passive  media  and  goes  beyond  the  eikonal  approximation 
of  geometric  optics  and/or  the  diffraction  theory  used  in  earlier  works  on  path  integrals 
in  electromagnetic  theory.  The  results  obtained  here  are  believed  to  be  useful  for  further 
development  of  path  integral  methods  in  theoretical  and  numerical  studies  of  propagation  of 
electromagnetic  wave  packets  in  passive  media. 

The  idea  of  numerical  simulations  follows  from  (1.7)  rather  straightforwardly,  namely, 

=  (1-8) 

Thus,  finding  a  state  of  the  system  in  a  sequential  moment  of  time  amounts  to  computing  the 
action  of  the  exponential  of  a  differential  operator  H  on  the  state  at  the  preceding  moment 
of  time.  Theoretically,  it  is  sufficient  to  know  ZYa*  up  to  (At)^.  The  limit  At  — >  0  in  (1.7) 
would  not  change  if  we  replace  the  exact  infinitesimal  evolution  operator  kernel  by  such  an 
approximation.  In  numerical  simulations,  the  limit  is  never  achieved.  Therefore  a  higher 
accuracy  is  required  to  make  errors  small.  Note  that  the  errors  are  accumulated  as  more 
iterations  (1.8)  are  taken.  An  expansion  of  exp(At7Y)  into  the  power  series  up  to  some 
desired  order  is  known  to  produce  unstable  algorithms.  Yet  another  obvious  drawback  is 
the  lack  of  imitaxity  of  the  time  evolution,  that  is,  if  the  Hamiltonian  is  skew-symmetric 
(anti-Hermitian,  if  a  complex  phase  space  is  used),  7Y*  =  —H,  then  U^Ut  =  1.  In  the 
Maxwell  theory,  as  we  shall  see,  the  squared  norm  of  with  respect  to  the  L2(]R®)  scalar 
product  (^’i.^’a)  =  / is  proportional  to  the  electromagnetic  energy  of  the  system. 
Consequently,  for  non-absorbing  media  the  unitarity  of  the  time  evolution  is  required  in 
simulations  to  provide  the  energy  conservation. 

We  shall  apply  Nelson’s  method  of  obtaining  the  path  integral  representation  of  the  fun¬ 
damental  solution  of  Maxwell’s  equations  for  passive  media.  It  is  based  on  the  Kato-Trotter 
product  formula  for  the  exponential  of  a  sum  of  two  noncommuting  operators  and  the  use 
of  the  Fourier  basis  to  compute  exponentials  of  differential  operators.  Actually,  in  practical 
applications,  a  simpler  version,  known  as  the  Lie-TVotter  product  formula,  is  used  (see  the 
textbooks  [13]  for  details  and  references  therein).  In  computational  quantum  mechanics  this 
is  also  known  as  the  split  operator  method.  It  allows  one  to  keep  the  differential  operators  in 
the  exponential,  and  thereby,  ensure  the  correct  dispersion  relation  of  simulated  electromag¬ 
netic  waves.  It  will  be  shown  that  there  exists  a  particular  realization  of  this  idea  in  which 
the  Gauss  law  holds  exactly  in  simulations.  In  general,  the  Gauss  law  can  be  enforced  by  the 
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projection  operator  formalism  developed  for  the  path  integral  representation  of  constrained 
d5mamical  systems  (for  a  review  see  [14]  and  references  therein,  a  numerical  application 
to  constrained  wave  packet  propagation  can  be  found  in  [15]).  The  idea  is  to  replace  the 
Hamiltonian  by  its  projection  on  the  subspace  (1.3).  If  P  is  the  projection  operator,  that 
is,  CaV^  =  0  for  any  —  T  and  V*  =  P,  then  H  is  replaced  by  VHV.  In  Maxwell’s 
theory  the  projection  can  be  implemented  in  our  algorithm  with  no  extra  computational 
costs.  A  significant  difference  from  the  quantum  mechanical  case  is  that  the  Hamiltonian 
W  (or  its  projection)  is  not  normal,  that  is,  it  does  not  commute  with  its  adjoint.  This 
feature  complicates  the  stability  analysis  because  the  von  Neumann  criteria  is  no  longer 
sufficient  for  stability,  while  still  being  necessary  [16].  Nevertheless,  the  stability,  accuracy 
and  convergence  analysis  of  the  algorithm  can  be  carried  out  in  rather  general  settings. 

Since  time  domain  simulations  are  performed  on  finite  lattices,  there  is  always  a  moment 
of  time  when  the  simulated  signal  first  reaches  the  lattice  boimdary.  One  typically  uses 
lattices  with  periodic  boundary  conditions.  So,  the  pulse  would  appear  on  the  other  side  of 
the  lattice  interfering  with  itself,  thus  leading  to  totally  disastrous  results  for  simulations. 
The  problem  is  usually  solved  by  introducing  absorbing  boundary  conditions  (see,  e.g.,  [17| 
(for  quantum  mechanics)  and  [18]  (for  electrodynamics)).  It  is  convenient  to  set  a  conducting 
layer  at  the  grid  boundary  whose  conductivity  is  chosen  so  that  it  neither  transmits  nor 
reflects  within  the  designated  accuracy  in  the  frequency  domain  of  the  initial  pulse.  In 
Appendix  we  briefly  describe  how  such  a  conducting  layer  can  be  obtained. 


2  Maxwell  theory  in  the  Hamiltonian  formalism 


Dynamics  of  electromagnetic  waves  in  continuous  media  is  governed  by  Maxwell’s  equations 

dtUt  =  cVxHt  (2.1) 

dtBt  =  -cVxEt,  (2.2) 

where  c  is  the  speed  of  light,  boldface  letters  denote  three-vector  fields  in  whose  spatial 
arguments  are  suppressed  and  the  time  dependence  is  indicated  by  a  subscript.  No  external 
currents  and  charges  (antennas)  are  included  in  this  study.  However,  the  formalism  being 
developed  is  readily  generalized  to  the  case  when  external  time  dependent  sources  are  present. 
The  electric  and  magnetic  induction  vectors,  Dt  and  Bf,  respectively,  are  subject  to  the 
constraints  (the  Gauss  law) 

V-Dt  =  V-Bt  =  0.  (2.3) 

In  linear  response  theory,  assumed  through  out  the  paper,  the  electric  induction  is  related 
to  the  electric  field  as  [19] 


Er  =  Et  -1-  Pj  , 


(2.4) 
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where  Xt  is  an  electric  response  function  of  the  medium  and  P*  is  the  medium  polarization 
vector.  A  similar  relation  can  be  written  for  the  magnetic  field  and  induction,  Bt  =  Ht+Mj, 
where  magnetization  Mj  is  determined  by  the  applied  magnetic  field  and  the  magnetic 
response  function  of  the  medimn. 

The  relation  between  inductions  and  fields  must  be  causal,  meaning  that  the  response  of 
the  medium,  P*  and  M(,  can  only  depend  on  fields  applied  to  the  medium  prior  to  the  current 
time  t,  (e.g.,  Xt  =  0  for  t  <  0)  [19].  A  natural  way  to  ensure  the  causality  is  to  require  that 
the  response  function  satisfies  a  differential  equation.  In  other  words,  the  response  function 
is  assumed  to  be  the  fundamental  solution  of  some  time  evolution  differential  equation. 
This  differential  equation  ean  be  obtained  from  a  particular  physical  model  of  the  medium 
in  question.  A  popular  model  is  the  multi-resonant  Lorentz  model.  Let  and  be 
the  Fourier  transforms  of  the  electric  induction  and  field.  Then  firom  (2.4)  it  follows  that 
D^,  =  The  dielectric  constant  in  the  Lorentz  model  has  the  form 


=  i+E 


pa 


,2  ^^2 


(2.6) 


and  Mf  =  0.  The  physical  meaning  of  the  model  is  transparent.  The  medium  is  assmned 
to  be  made  of  N  sorts  of  damped  harmonic  oscillators  with  resonant  fi:equencies  uJa  and 
damping  coefficients  7o.  Parameters  Upa,  called  the  plasma  frequencies,  are  proportional  to 
coupling  constants  of  the  oscillators  to  the  external  electric  field  (the  electric  dipole  coupling) 
and  also  depend  on  the  density  of  oscillators  of  the  sort  a.  The  density  may  vary  in  space. 
So  u)pa  are  assumed  to  be  functions  of  spatial  coordinates.  In  an  empty  space,  Upa  =  0.  If 
the  resonant  frequency  is  zero,  the  one-resonant  Lorentz  model  is  equivalent  to  the  Drude 
model  of  metals  [19].  In  the  Lorentz  medium  the  magnetic  response  function  is  zero,  while 
the  electric  response  function  can  easily  be  found  by  taking  the  Fourier  transform  of  (2.5). 
Its  expheit  form  is  omitted  here  because  it  will  not  be  used.  The  medimn  polarization  is 
determined  by  a  set  of  second-order  differential  equations 


N 


-h  27«atP?  -h  ulPt  =  ,  P*  =  P? 


(2.6) 


a=l 


Together  with  Maxwell’s  equations,  Eq.  (2.6)  form  a  system  of  sought-for  causal  evolution 
equations  which  are  to  be  transformed  into  a  system  of  first  order  equations  by  means 
of  the  Hamiltonian  formalism.  In  finite  difference  time  domain  numerical  schemes,  the 
Hamiltonian  formalism  for  the  Lorentz  model  has  been  used  in  [20]  to  study  propagation  of 
an  electromagnetic  pulse  in  homogeneous  Lorentz  media. 

In  our  approach  no  boundary  conditions  arc  imposed  on  electromagnetic  fields  at  medium 
and/or  target  interfaces.  The  latter  are  modeled  by  spatially  dependent  couplings  of  media  to 
electromagnetic  fields  which  are  included  into  the  system  Hamiltonian.  At  any  interface,  the 
couplings  are  allowed  to  have  discontinuities,  or,  from  a  physical  point  of  view,  they  remain 
smooth  but  change  rapidly,  \w\VijJp\lu)p  :s>  1,  at  the  interface,  where  A^,  is  a  typical  wave 
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length  of  the  incoming  wave  packet.  The  conventional  boundary  conditions  are  automatically 
generated  by  the  dynamics  [19].  Thus,  the  initial  value  problem  is  solved  in  L2(1R^)  for  every 
matter  and  electromagnetic  field  component.  This  implies  that  the  energy  of  the  propagating 
wave  packet  remains  finite  (in  contrast  to  the  scattering  matrix  approach  based  on  plane 
wave  solutions). 

Let  us  now  formulate  the  initial  value  problem  for  a  generic  passive  medium  and  then 
apply  the  formalism  to  multi-resonant  Lorentz  models.  Combine  the  fields,  inductions  and 
medium  responses  into  columns: 


Assuming  linear  response  theory,  one  can  write  for  the  Fourier  transforms 

5  (2-7) 

where  the  Fourier  transform  of  a  general  response  function,  Xa>,  has  to  satisfy  a  dispersion 
relation  that  ensmres  causality  (like  the  Kramer-Kronig  relations  for  the  dielectric  constant) 
[19].  For  anisotropic  media,  Xw  is  a  symmetric  matrix  acting  on  components  of  electromag¬ 
netic  fields.  With  this  type  of  generality  all  possible  media  are  covered  as  long  as  linear 
response  theory  is  valid.  The  response  function  Xw  can  either  be  modeled  or  measured  and 
tabulated  in  some  firequency  range  of  interest  (determined  by  the  frequency  bandwidth  of 
the  initial  wavepacket),  say,  oj  €  [a;i,a;2].  Next,  the  components  of  expanded  in 

a  basis  of  suitable  orthogonal  polynomials.  An  optimal  expansion  is  often  achieved  in  t^e 
Chebyshev  polynomial  basis.  Chebyshev  polynomials  are  defined  in  the  interval  [—1, 1]  so 
a  corresponding  rescaling  and  translation  of  [wi,a;2]  must  be  done.  By  taking  the  Fourier 
transform  of  xZ^'^w  ~  obtain  the  desired  differential  equation 

N 

JZxnSr^f  =  WpV’f  ,  (2.8) 

n=0 

where  Up  =  a;p(r)  plays  the  role  of  the  coupling  constant  between  matter  and  electromagnetic 
fields.  The  order  N  is  determined  by  the  highest  order  of  polynomials  used  to  approximate 
The  expansion  coefficients  Xn  and  the  couphng  Wp  are  matric^  for  anisotropic  media. 

The  basic  idea  of  the  Hamiltonian  formaUsm  is  to  convert  the  system  (2.6)  or  (2.8)  into 
a  system  of  first-order  differential  equations  by  introducing  auxiliary  (matter)  fields.  The 
nmnber  of  such  fields  is  determined  by  the  order  of  the  original  evolution  equation  for  matter. 
For  instance,  in  the  case  of  the  multi-resonant  Lorentz  model,  there  are  N  fields  P“  each  of 
which  satisfies  a  second  order  differential  equation.  In  the  Hamiltonian  formalism  one  would 
have  2N  real  vector  fields,  j  =  1, 2, ...,  2N.  A  simple  possibility  is  to  set 


P“ 

(2.9) 

= 

(2.10) 

=  ~27a^t'*  —  °  ^  +  ^pa^t  • 

(2.11) 
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The  reasor  of  inserting  the  factor  uipa/^a  fte  definition  (2.9)  of  the  auxihary  fields  will 
be  evident  from  what  follows.  Note  that  the  medium  polarization  P*  must  be  zero  in  empty 
space  where  Upa  =  0.  The  factor  in  (2.9)  simplifies  the  energy  conservation  and  stability 
analysis. 

For  the  Lorentz  model  there  is  another  convenient  way  to  introduce  the  Hamiltonian  for¬ 
malism  by  using  N  complex  vector  fields  which  satisfy  the  first  order  differential  equation 

dtCt  =  (2.12) 

p?  =  ^  (c + ca .  p-13) 

where  =  -7a  +  iva  and  Va  =  -  72.  This  representation  is  defined  only  if  7a  <  Wa 

(i.e.,  the  attenuation  is  not  high).  Prom  the  numerical  point  of  view,  solving  a  decoupled 
system  of  N  first  order  differential  equation  and  taking  complex  conjugation  (denoted  here 
by  an  over  bar)  is  less  expensive  than  solving  an  original  system  of  differential  equations  to 
compute  the  medium  polarization. 

Returning  to  the  general  case,  we  introduce  a  set  of  auxihary  fields  to  convert  (2.8) 
into  a  first-order  system, 

dt^t^nUt  +  VMF^[  .  (2.14) 

The  operators  and  Vmf  are  determined  by  the  details  of  going  over  to  the  Hamiltonian 
formalism.  We  shall  call  Hm  the  matter  Hamiltonian;  it  governs  time  evolution  of  the 
medium  when  no  external  fields  are  applied.  The  index  F  indicates  that  the  electromagnetic 
degrees  of  freedom  are  described  by  fields,  not  inductions.  We  shall  see  shortly  that  the 
matter  Hamiltonian  depends  on  whether  or  ipf  is  used  as  independent  electromagnetic 
variables.  The  matrix  Vmf  describes  the  coupfing  of  matter  to  the  electromagnetic  fields, 
which  is  emphasized  by  the  index  MF  (matter-to-field  coupfing).  We  introduce  a  linear  time 
independent  operator  TZ  that  acts  in  the  space  of  auxiliary  (matter)  fields  so  that 

=  (2.15) 

that  is,  the  (response)  operator  1Z  maps  a  given  configmation  of  auxiliary  fields  onto  the 
corresponding  physical  response  field.  It  depends  on  the  definition  of  the  matter  fields  (cf. 
(2.9)  and  (2.13)).  A  passive  medium  is  not  excited,  =  0,  if  no  external  electromagnetic 
field  is  applied;  that  is,  the  initial  condition  for  Eq.  (2.14)  is  such  that  it  has  only  the  trivial 
solution  whenever  ipf  —  0.  Under  this  condition,  the  solution  of  (2.14)  reads 

J —00 

Hence,  the  linear  response  operator  Xu}  in  (2.7)  is  the  Fourier  transform  of  the  operator 
Xt  =  6t'Rexp{tHM)VMF,  where  6t  is  the  Heaviside  function.  Or,  vice  versa,  11,  Hm  and 
Vmf  must  be  chosen  so  that  the  Fourier  transform  of  Xt  coincides  with  the  known  response 
function  Xa>  of  the  medium  in  a  designated  frequency  range. 
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Maxwell’s  equations  without  external  currents  can  be  rewritten  in  the  Hanultonian  form 

dtrl>f  =  HFi>f-dtipf  =  HFi’f  +  VFMit.  (2.16) 

The  field-tomatter  coupling  Vfm  and  the  field  Hamiltonian  Hf  are  deduced  firom  (2.14)  by- 
acting  on  the  latter  by  the  operator  'R.,  which  jdelds 

Vfm  =  >  (2-17) 

Hf  =  ^  “  TJVmf  =  7^0  “  ^Vmf  •  (2-18) 

It  is  always  possible  to  set  up  the  Hanultonian  formalism  so  that  RVmf  —  0  and,  hence, 
Hf  =  Ho-  It  is  not  difficult  to  verify  that  this  holds  for  the  Lorentz  model  discussed  above. 
In  the  general  case,  the  standard  procedure  of  going  over  to  the  Hamiltonian  formalism  [8], 
where  components  of  are  identified  with  time  derivatives  of  the  response  field,  ~  > 

leads  to  the  same  result  that  HVmf  =  0.  Thus,  without  loss  of  generality,  the  last  term  in 
the  field  Hamiltonian  (2.18)  can  be  omitted. 


The  auxiliary  matter  and  electromagnetic  fields  (or  inductions)  are  unified  into  a  larger 
coliunn  ,  „ , 


The  wave  function  satisfies  the  Schrodinger  equation 


)■ 

(2.19) 

Vfm^ 

nLJ  ■ 

(2.20) 

which  has  to  be  solved  with  the  initial  field  configuration  while  the  matter  fields 

are  assumed  to  be  zero  at  the  initial  moment  of  time,  ^t=o  =  0,  e.g.,  the  initial  wave  packet 
is  localized  in  an  empty  space  region.  Equations  (2.14)  and  (2.16)  are  equivalent  to  (2.20). 
In  a  similar  fashion,  one  can  derive  the  Schrodinger  equation  for  .  Note  that 


5=(j  T)'  ^'=(j  "r)> 

Hence, 

dti^l  =  H^^l  ,  H^  =  SH^S-^  .  (2.22) 

The  corresponding  blocks  of  H^  have  the  form 

Hi  =  Ho  ,  Vmi  =  Vmf  >  (2.23) 

Vm  =  Vfm +  nHM-Hon  =  -Hon,  (2-24) 

Hii  -  Hfi-VMFn.  (2.25) 


To  simplify  Vm,  Eq.  (2.17)  has  been  applied.  Observe  in  (2.25)  the  aforementioned  de¬ 
pendence  of  the  matter  Hamiltonian  on  the  representation  of  electromagnetic  degrees  of 
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freedom.  The  use  of  either  (2.20)  or  (2.22)  in  numerical  simulations  has  its  own  advantages 
and  disadvantages  which  are  discussed  below. 

As  an  example,  we  give  an  explicit  form  of  the  Hamiltonian  for  the  Lorentz  model  when 
the  auxiliary  field  are  defined  by  (2.9) 

(2.26) 

(2.27) 

(2.28) 

where  diag  indicates  that  the  corresponding  matrix  is  block-diagonal  with  blocks  listed  in 
the  order  from  the  upper  left  to  lower  right  corners.  Note  that  the  matrices  VpMa  aJid 
UpMa  act  on  a  ^-dimensional  column  Therefore  they  should  be  understood 

as  composed  of  3  x  3  blocks.  Each  block  is  obtained  by  multiplying  the  unit  matrix  by  the 
munber  indicated  in  place  of  the  block  in  (2.26)  and  (2.28). 

Our  final  remark  in  this  section  concerns  “canonical”  transformations  in  the  Hamiltonian 
formalism.  As  has  been  pointed  out,  the  auxiliary  fields  arc  not  uniquely  defined.  There  is 
a  freedom  of  making  general  complex  nonsingular  linear  transformations  such  as 

det5M?^0.  (2.29) 

If  the  infinitesimal  evolution  operator  =  exp(Ai7f(^’^))  is  computed  with  one  choice  of 
the  auxiliary  fields,  a  simple  similarity  transformation,  like  the  one  in  (2.22),  would  allow  us 
to  compute  it  in  any  other  basis  of  auxiliary  fields.  This  is  an  important  observation  because 
the  auxiliary  field  basis  can  be  chosen  in  a  way  that  facilitates  computation  of  the  evolution 
operator  (e.g.,  to  improve  the  convergence  rate  or  speed  up  simulations).  For  instance,  in 
the  complex  representation  (2.13)  of  the  auxiliary  fields  in  the  Lorentz  model,  the  matter 
Hamiltonian  is  diagonal.  The  corresponding  transformation  of  the  auxiliary  fields  is  given 

by 


To  transform  the  whole  system  into  this  representation,  the  Hamiltonian  is  replaced  by 
S~^H^S  and  the  wave  function  by  where  S  is  block-diagonal  with  the  unit  matrix 
in  the  upper  left  (field)  comer  and  with  Sm  in  the  lower  right  (matter)  corner. 


Vfm 

Va/f 

'ijF 


=  {Vfm 
=  -VI 


1,VfM2,  '  •  •  ,VfMn)  ,  VFMa=^Q  j 


—  ~  ) 

=  diag  {Hmv  i'^Mn)  >  —27a)  ’ 


3  The  grid  representation  of  Maxwell’s  theory 

Consider  an  equidistantly  spaced  finite  grid  with  periodic  boundary  conditions.  Let  Ar 
be  the  grid  step  and  n  be  a  vector  with  integer  valued  components.  Then  the  dynamical 
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variables  are  projected  onto  the  grid  by  taking  their  values  at  grid  points  r  =  nAr , 

’J'f  (r)  li'P(nAr)  ,  (3.1) 

where  Q  denotes  the  representation,  I  or  F.  For  simplicity,  a  cubic  grid  is  assumed  here. 
It  is  straightforward  to  generalize  the  discussion  to  a  generic  rectangulax  grid.  Consider  a 
discrete  Fourier  transformation  associated  with  the  grid  [5,  21] 

|rf(nifco)  =  X;-^nn''J'?(n'Ar)  ,  F*F  =  Fjr  =  l,  (3.2) 

n' 

where  the  dual  lattice  step  is  fco  =  27r/Ar.  The  grid  spatial  size  L  and  step  must  be 
chosen  so  that  the  Fourier  transform  of  the  initial  wavepacket  has  support  within  the  region 
k  e  [kmin,kmax]  where  A;  =  |k|,  k„u^  =  k)  and  kmin  =  The  Hamatonian  'HP  is  split 

into  a  sum 

+  (3.3) 

where  all  the  spatial  derivatives  are  included  into  ’Hq  ^  contains  multiplications  by 
position  dependent  functions.  This  is  always  possible  for  the  Hanultonian  described  in  the 
preceding  section.  The  operator  V®  is  projected  naturaUy 

V^(r)^f  (r)  -+  V^(nAr)^ff  (nAr)  .  (3.4) 

Consider  in  the  Fourier  basis,  W^(V)  ->  W^(ik).  The  projection  is  then  done  via  the 
discrete  Fourier  transform 

Wo(V)’^f(r)|  (3.5) 

!r±=nAr 

n' 

The  projection  (3.5)  as  well  as  any  action  of  on  state  vectors  is  performed  by  the  fast 
Fourier  method  [5].  It  requires  IVloga  N  elementary  operations  (flops)  with  N  being  the  grid 
size. 

In  what  follows,  the  rules  (3.4)  and  (3.5)  define  the  action  of  the  operators  and  Hq 
and  their  functions  on  any  state  vector.  The  action  of  a  product  of  and  on  any  state 
vector  is  imderstood  as  consecutive  actions  of  these  operators  according  to  the  rules  (3.4) 
and  (3.5),  in  the  order  specified  in  the  product. 

In  finite  differencing  schemes,  the  action  of  Hq  on  a  state  vector  would  require  mNd  op¬ 
erations  where  the  integer  m  depends  on  a  particular  difference  scheme  used  to  approximate 
derivatives,  and  Nd  is  the  grid  size  used  in  the  differencing  scheme.  It  should  be  noted  that, 
as  shown  below,  the  use  of  the  fast  Fourier  transform  eliminates  the  phase  error  (because  the 
correct  electromagnetic  dispersion  relation  is  preserved)  and  operates  at  the  Nyquist  limit. 
These  two  features  allows  one  to  reduce  substantially  the  grid  size  as  compared  with  that  in  a 
finite  differencing  scheme,  while  providing  the  same  accuracy  in  simulations.  Recall  that,  in 
scattering  problems,  the  phase  of  the  return  signal  contains  the  most  significant  information 
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about  the  target.  So,  in  practice,  grids  in  finite  differencing  schemes  are  significantly  larger 
(more  dense)  than  grids  used  in  the  fast  Fourier  method  in  order  to  reduce  the  phase  errors 
in  the  former.  Needless  to  say,  the  advantage  of  the  fast  Fourier  method  in  reducing  the 
phase  error  becomes  even  more  significant  in  higher  dimensions  because  Nd/N  =  {ud/n)^ 
where  rid  and  n  are  the  corresponding  numbers  of  grid  points  per  shortest  wave  length  in 
the  wave  packet,  and  D  is  the  grid  dunension.  The  Nyquist  limit  is  n  =  2,  while  nd  is  of 
order  10  or  higher. 


4  The  split  operator  method 

Let  I'i'l  denote  the  L2(M®)  norm  of  the  wave  function,  or  the  Eucfidean  norm  of  the  corre¬ 
sponding  vector  (3.1)  in  the  grid  representation.  One  possible  way  to  compute  nxunerically 
the  path  integral  (1.7)  is  based  on  the  Kato-Iirotter  product  formula  [13] 

■  lim  e*(-4+S)^_  (eW2ng*S/ngt4/2n)"^  (4.1) 

n~+oo 

for  a  general  ^  and  under  certain  assumptions  about  the  linear  operators  and  B  in  the 
Hilbert  space  spanned  by  For  our  purposes  it  is  sufficient  to  note  that  for  bounded 
operators,  (4.1)  always  holds  and  is  known  as  the  Lie- Trotter  product  formula.  In  the  grid 
representation,  which  would  always  be  assumed,  unless  stated  otherwise,  operators  A  and  B 
are  finite  matrices  and,  hence,  bounded. 

Let  us  apply  (4.1)  to  the  split  (3.3),  meaning  that  the  operator  Hq  is  used  in  plane  of 
A  (or  B)  and,  respectively,  the  operator  is  used  in  place  of  B  (or  A).  The  infinitesimal 
evolution  operator  in  (1.8)  can  be  approximated  by  the  first  term  in  the  following  expansion 


uS, 

(4.2) 

rQ 

y^t 

^  ^AtA/2^AtB^AtA/2 

(4.3) 

=  -  ^  (lA  lA  ell  -  2[B,  [B.  .4)1)  +  0(At)  . 

(4.4) 

By  making  n  larger  while  keeping  nAt  =  t  fixed,  the  strong  convergence  in  (4.1)  guarantees 
that  the  error  can  be  made  arbitrary  small  for  any  initial  state, 

|(ZYf,)"^o-(a^,)”'^o|--0  (4.5) 

as  n  — oo  for  any  'i'o-  The  numerical  iteration  algorithm  is  then  based  on  the  replacement 
of  the  exact  evolution  (1.8)  by  the  approximate  one 

(«) 

The  quantity  tAt®|VV^$ol/l®ol  can  be  used  to  roughly  estimate  the  accuracy  of  the  al¬ 
gorithm.  A  more  detailed  accuracy  analysis  is  given  in  Section  8.  By  making  use  of  the 
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Campbell-Hausdorf  formula  for  the  exponential  of  the  sum  of  operators  it  is  possible  to  ob¬ 
tain  the  symmetric  product  formula  in  (4.2)  to  approximate  lAt^i  up  to  any  desired  order 
in  At,  presumably  achieving  a  higher  accuracy  [22].  This  would  come  at  the  price  of  hav¬ 
ing  more  exponentials  in  the  symmetric  product  In  numerical  simulations,  one  should 
keep  in  mind  that  computational  costs  of  decreasing  At  in  the  third  order  split  (4.2)  (i.e., 
increasing  the  number  of  steps  in  the  time  evolution)  might  be  less  than  those  of  computing 
a  lesser  number  of  actions  of  in  higher  order  splits.  So,  the  higher  order  splits  are  not 
always  optimal  to  achieve  a  better  accuracy  [2] . 


On  the  grid,  the  action  of  the  amplification  operator  Q%.  is  computed  according  to  the 
rules  (3.4)  and  (3.5)  applied  to,  respectively,  cxp(AtV'^)  and  cxp(Am^).  ExpUcit  formulas 
for  the  exponentials  of  the  corresponding  operators  can  be  worked  out  in  the  field  and 
induction  representations.  If  the  fields  are  used  as  independent  variables,  then  a  natural 
choice  is 


^.)  +  ( 


0 

.  Vmf 


Vfm\ 

0  ) 


(4.7) 


The  matter  Hamatonian  Hi,  can  also  be  transferred  into  if  so  desired.  This  rearrange¬ 
ment  affects  the  accuracy  of  the  method,  meaning  that  the  operator  (4.4)  would  change.  In 
turn,  a  rearrangement  of  operators  in  the  split  can  be  used  to  improve  the  accuracy.  We 
shall  discussed  this  issue  later.  Using  the  Taylor  series  we  infer  that 


exp(fWo)  =  1  +  ^cos(ctv'-A)  -  l]  Pj.  +  c\/^A~ 

where  A  =  V  •  V  is  the  Laplace  operator  and  'Px  =  1  —  V(A~^V'  )  is  the  projector 
on  transverse  fields,  that  is,  P±E  =  E  if  V  •  E  =  0  and  Pj^E  =  0  if  the  vector  field 
E  is  conservative,  E  =  V^.  The  projector  P±  can  be  omitted  in  (4.9)  if  it  is  known 
(e.g.  from  a  theoretical  analysis  of  the  system)  that  the  fields  remain  transversal  in  due 
course.  In  this  case,  the  two  first  terms  in  (4.9)  are  equal  to  cos(ct\/— A).  The  action  of 
exp(tHo)  is  computed  by  the  fast  Fourier  transform  according  to  (3.5).  In  the  Fourier  basis, 
—A  =  =  n^kl-  Note  also  that  the  Fourier  transform  of  the  fields  is  required,  while 

the  auxiliary  fields  remain  in  the  grid  basis  all  the  time.  The  exponentials  of  Hi,  and 
can  either  be  computed  analytically  for  simple  models  like  the  Lorentz  model,  as  is  shown 
Section  5,  or,  in  general  case,  by  direct  diagonaUzation  at  each  grid  site. 


Alternatively,  the  following  approximation  of  the  exponential  of  an  operator  can  be  used 


1  -h  AtB/2  +  At^gVl2  5 

1  -  AtB/2  +  At2BVl2  ^  ^ 


1  +  AtB/2 
1  -  AtB/2 


+  0{At^) . 


(4.10) 


If  the  matrix  B  is  anti-hermitian,  then  the  approximations  (4.10)  of  the  exponential  of  B 
retain  unitarity,  which  is  important  for  stability  of  the  split  algorithm  (see  Theorems  7.1 
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and  7.3).  On  the  other  hand,  costs  of  computing  the  inverse  matrices  in  the  right  hand  side 
of  (4.10)  can  be  less  than  those  of  computing  the  exponential. 


If  the  inductions  are  used  as  independent  variables,  then  a  natural  choice  of  the  split 
would  be  /  X  /  \ 

4^).  (4.11) 


Making  use  of  the  Taylor  expansion  again  we  deduce  that 


J 


and,  similarly. 


exp(fV  )  -  1]Vmj  expitUif) 


Af)) 


(4.12) 


(4.13) 


Now  we  can  compare  the  two  splits.  The  split  (4.7)  has  an  advantage  over  (4.11)  because 
it  requires  less  calls  of  the  fast  Fourier  transform.  Indeed,  in  the  former  the  fast  Fourier 
transform  is  called  only  for  the  fields  V’f  •  As  it  follows  firom  (4.12),  the  operator  exp(t'Ho) 
acts  on  both  the  inductions  V’f  and  the  auxiUary  fields.  Hence  the  fast  Fourier  transform  must 
be  called  for  the  entire  column  .  If  the  number  of  auxihary  fields  is  large,  there  might  be 
a  substantial  difference  in  the  computational  speed  of  two  algorithms.  The  latter,  however, 
depends  on  the  choice  of  the  matter  fields  which,  in  turn,  determines  TZ.  and  therefore  the 
number  of  flails  of  the  fast  Fourier  transform.  Note  that  the  matter  fields  can  always  be 
chosen  in  such  a  way  that  only  one  of  the  components  of  specifics  the  response  field  . 
Thus,  the  canonical  transformation  (2.29)  can  be  used  to  reduce  the  number  of  calls  of  the 
fast  Fourier  transform.  If  H  is  chosen  so  that  it  depends  on  position,  multiplication  of  the 
matter  fields  by  TL  must  be  done  before  calling  the  fast  Fourier  transform.  A  significant 
advantage  of  the  split  in  the  induction  representation  is  that  the  Gauss  law  can  be  exactly 
fulfilled  without  altering  the  algorithm  (see  Theorem  8.2). 

In  empty  space  either  of  the  splits  reproduces  an  exact  solution  of  Maxwell’s  equations 
for  any  period  of  time  t,  provided  the  initial  pulse  is  bandwidth  limited.  Indeed,  on  the 
grid,  the  initial  wave  packet  is  a  superposition  of  a  finite  number  of  plane  waves.  Thanks  to 
the  linearity  of  the  theory,  each  Fourier  mode  is  evolved  exactly,  that  is,  without  any  phase 
error,  by  exp{tHo)  for  any  t  >  0.  As  final  remarks  in  this  section,  we  note  that  the  algorithm 
can  operate  at  the  Nyquist  fimit:  Two  grid  points  per  shortest  wavelength  in  the  initial 
wave  packet  [5].  Yet,  for  the  multiresonant  Lorentz  model,  it  is  unconditionally  stable  (see 
Theorems  7.1  and  7.3).  These  features  cannot  be  achieved  in  any  finite  difference  scheme. 
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5  A  multi-resonant  Lorentz  model 

An  analytical  expression  for  the  exp>onents  of  and  helps  to  reduce  computational 
costs.  Here  such  anal3rtical  expressions  are  derived  for  multi-resonant  Lorentz  models.  Let 
us  take  first  the  field  representation.  Due  to  the  block  diagonal  structure  of  we  get 

exp(tW^)  =  diag  (exp(tW^i),  exp(t7i^2)>  ”*»  exp(t?<^jv))  ,  (5-14) 

exp(iW^„)  =  e"'*'”*  cosh  Uat  +  ^  +  7a)  »  (5.15) 

where  Pa  =  ill  -  The  exponential  (5.15)  is  easy  to  compute  by  expanding  in 

the  Pauli  matrix  basis,  which  is  also  a  basis  for  the  Lie  algebra  su(2),  and  then  by  using 
the  well  known  formula  for  the  exponential  of  a  hnear  combination  of  Pauli  matrices.  For 
small  attenuation,  7a  <  Wa,  we  get  Pa  =  iva-  The  hyperbolic  fimctions  in  (5.15)  become 
trigonometric  ones  and  Pa  is  replaced  by  !/„•  The  eigenvalues  of  the  matter  Hamiltonian  are 
Aa  =  — 7a±i?a.  Hence,  Re  Aa  <  0  and  amplitudes  of  the  matter  fields  are  always  exponentially 
attenuated  as  t  — >  oo,  unless  7a  =  0  leading  to  Re  Aa  =  0. 

Computation  of  exp(tV^)  is  a  bit  more  subtle.  We  start  with  the  observation  that  the 
characteristic  polynomial  of  has  a  simple  form 

det  (V'’  -  A)  =  A’-'CA^  +c4).  = 

0=1 

This  can  be  proved  either  by  a  direct  computation  or  by  mathematical  induction.  So, 
has  2N  zero  eigenvalues  and  two  non-zero  ones,  A  =  itWp.  Let  X  be  the  eigenvector  of 
corresponding  to  the  eigenvalue  iujp.  Its  components  have  the  form 

=  s'  [Hi  +  J  2, 2(iV -t- 1)  , 

so  that  X  ‘  X  =  1  and  X  •  X  =  X  •  X  =  0  where  the  dot  denotes  the  Euclidean  scalar 
product.  The  skew-symmetric  matrix  has  the  following  spectral  decomposition 

=^iu}p{X®X -X0X)  .  (5.17) 

Taking  the  square  of  (5.17)  we  also  infer  that 

X  ®  X  = .  (5.18) 

The  exponential  of  (5.17)  is  obtained  via  the  Taylor  series  and  making  use  of  (5.18).  The 
final  result  reads 

^(tVn  =  l  +  5^V''  +  2f2!(!VZ2)y  .  (5.19) 

S  \  S  / 
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In  the  induction  representation,  an  explicit  formula  for  exp(t'Hj[^)  is  not  that  simple. 
To  avoid  unnecessary  technicalities,  we  limit  the  discussion  to  the  simplest  case  of  the  on^ 
resonant  Lorentz  model.  We  choose  the  matter  fields  so  that  =  Pt  ~  ^  this 

case,  non-zero  elements  of  the  matter  Hamiltonian  are  ~  ~  ~  ^1’ 

Hif22  =  -27.  The  coupling  mattrix  Vmi  has  only  one  non-zero  element,  Vum  =  Here 
0^0  is  the  resonant  frequency,  7  is  the  attenuation  constant  and  Wp  is  the  plasma  frequency. 
Using  the  Pauli  matrix  basis  again,  we  find  that  the  expression  (5.15)  holds  for  expitHif)  if 
we  replace  in  it  i/„  by  Vp  =  y'wg -1- -  72,  7„  by  7  and  Hf,  by  The  lower  left  comer 
of  (4.13)  has  the  form 

(Wl,)"'(exp(t?fif)-l)VM/  =  -;;|^(exp(<Wi,)-l)  (j  q)  •  (5-20) 

A  further  simplification  can  be  achieved  by  going  over  to  the  complex  representation  (2.12) 
of  the  auxihary  fields  in  which  the  matter  Hamiltonian  is  diagonal.  The  transformation  rule 
is  explained  in  the  paragraph  after  Eq.  (2.30). 


6  Energy  and  norm  conservation 

Consider  the  1L2(E®)  norm  of  ).  In  the  grid  representar 

tion,  the  norm  coincides  with  the  corresponding  (complrac)  Euclidean  norm,  ,  ^^)  =  Ar 
Yin  (^^  Ar )  vl/^(n  Ar)  where  the  sum  is  taken  over  all  grid  sites.  By  taking  the  time  deriva^ 

tive  of  P  and  using  the  evolution  equation  (1.1),  it  is  not  hard  to  deduce  that  the  norm 
is  conserved,  provided  the  Hamiltonian  is  anti-Hermitian 

=  -n^ .  (6.1)  ■ 

For  a  generic  passive  media  this  is  not  the  case.  So  the  norm  is  generally  not  conserved 
in  contrast  to  the  quantum  mechanical  case.  However,  we  shall  see  that  in  the  case  when 
the  matter  evolution  is  described  by  second  order  differential  equations  in  time  and  no 
attenuation  is  present,  the  norm  coincides  with  the  system  energy  and  is  conserved.  In 
numerical  simulations,  this  important  property  can  be  used  to  help  to  control  the  accmacy. 

Consider  multi-resonant  Lorentz  models  with  no  attenuation,  70  =  0.  We  start  with 
the  observation  that  the  field  and  matter  evolution  equations  can  be  obtained  from  the 
variational  principle  for  the  action 

,  (6.2) 

where  the  electromagnetic  degrees  of  freedom  are  described  by  vector  and  scalar  potentials, 
respectively.  At  and  <pt,  so  that  Et  =  Vtpt  -  and  Bf  V  x  At.  The  units  are  chosen  in 


S  =  y  i  (B?  -  B?)  +  ^  5;  ■  E, 
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this  Section  so  that  c—l.  The  polarization  of  the  medium  is  expressed  via  the  matter  fields 
as  Pt  =  Yla^pa^t-  The  least  action  principle  for  the  scalar  potential  (pt  leads  to  the  Gauss 
law,  V  •  Dt  =  0,  for  the  vector  potential  At  to  the  Maxwell’s  equation,  dtDt  =  V  x  Bt, 
and  for  the  matter  field  to  the  medixun  polarization  evolution  equation  of  the  Lorentz 
model  with  no  attenuation,  7o  =  0.  The  second  Maxwell’s  equation  and  the  Gauss  law  for 
the  magnetic  field  follow  from  the  relation  B*  =  V  x  At  by  taking  its  time  derivative  and 
divergence,  respectively.  The  energy  of  the  system  coincides  with  the  canonical  Hamiltonian 
which  is  obtained  by  a  Legendre  transformation  [8]  of  the  Lagrangian  L  for  the  velocities  ^  A* 
and  Doing  the  Legendre  transformation,  we  find  the  canonical  Hamiltonian  (energy) 

of  the  systena 

£,  =  i  I*  [e?  +  B?  +  ;^  «  +  ul^)  =  i  |«f  r  .  (6.3) 

where  w®  =  SL/S(dt'd^)  arc  canonical  momenta  of  the  matter  fields.  To  get  the  last  equality 
in  (6.3),  we  have  used  the  relations  and  which  follow  from  comparison 

of  the  canonical  Hamiltonian  equations  of  motion  for  the  canonically  conjugate  variables 
and  TT?  and  Eqs.  (2.10)  and  (2.11)  with  7„  =  0.  Note  that  the  canonical  momentum 
conjugate  to  the  vector  potential  At  coincides  with  — Dt  =  — Et  —  Pt,  not  — Et  in  this 
system.  Therefore,  the  coupling  between  the  electromagnetic  and  matter  degrees  of  freedom 
is  included  into  the  term  =  (Dt  -  Pt)^  of  the  canonical  Hamiltonian.  Equation  (6.3) 
becomes  the  conventional  expression  for  the  electromagnetic  energy  in  a  passive  medium 
[19]  when  irf  and  are  replaced  by  the  corresponding  solutions  of  the  equations  of  motion 
with  initial  conditions  Trg  =  =  0-  The  energy  conservation  can  be  deduced  either  from 

the  Noether  theorem  (because  Et  is  the  Noether  integral  of  motion  corresponding  to  the 
time  translational  symmetry  of  the  action)  or  directly  from  the  norm  conservation  of 
(because  the  evolution  operator  exp{tH^)  is  unitary  when  7a  =  0). 

In  numerical  simulations,  an  exact  unitary  evolution  operator  U^t  is  replaced  by  its 
approximation  However,  the  energy  remains  conservative: 

Theorem  6.1.  The  split  algorithm  is  imitary  for  multiresonant  Lorentz  models  with  no 
attenuation,  that  is,  the  spUt  algorithm  preserves  the  energy  Et+At  =  Et- 

Proof.  In  the  field  representation.  Ho*  =  -Hq  and  V^*  =  -V^  and,  therefore,  Gm  is 
unitary.  As  a  result,  the  algorithm  preserves  the  initial  wave  packet  energy  and  the  norm, 

I  =  Ka,I  =  K I .  {6-4) 

In  the  induction  representation,  the  energy  coincides  with  the  norm  of  in  the  measure 
space.  The  measure  is  determined  by  the  transformation  law  , 

=  =  =  (6.5) 

Since  H^  is  similar  to  H^,  the  Hamiltonian  H^  is  anti-Hermitian  relative  to  the  fi  scalar 
product, 

=  -^iH^  .  (6.6) 
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The  norm  conservation  (imitarity)  in  the  split  algorithm  requires  in  addition  that  the  am¬ 
plification  matrix  ei,  satisfies  the  imitarity  condition 

=  (6.7) 

This  is  indeed  the  case.  To  prove  (6.7),  we  show  that  Hq  and  satisfy  the  condition  (6.6) 
and,  hence,  the  product  of  their  exponentials  is  a  unitary  operator  relative  to  the  n  scalar 
product.  Consider  =  S'H^S~^  =  Til  +  “  S~^TilS  -4-  S~^V^S.  For  the 

Lorentz  model, 

5-‘74S=(’J°  o)=-(‘S'‘WJ‘S)*  .  (6.8) 

Therefore  Til  satisfies  (6.6).  Prom  the  anti-Hermiticity  of  Ti^  and  (6.8)  it  follows  that 

{S-^V^Sy  =  -S-^V^S  .  (6.9) 

Hence,  also  satisfies  (6.6).  Thus, 

= i»;i„ .  (6.10) 

which  completes  the  proof. 

The  norm  (energy)  conservation  can  be  used  to  control  numerical  convergence,  especially 
when  the  aliasing  problem  in  the  fast  Fourier  transform  is  present,  i.e.,  when  parameters 
of  the  medium  are  discontinuous  functions  in  space.  In  a  properly  designed  algorithm  the 
loss  of  energy  (norm)  due  to  attenuation  should  be  controlled  by  the  symmetric  part  of  the 
Hamiltonian  operator 

dtEt  =  I"*  =  \  (a-?,  <  0  ,  (6.11) 

a 

where  =  V^  =  {Ti^*  -f  'H9)/2  <  0  (a  negative  scmidcfinitc  operator)  which  is,  in  this 
case,  a  diagonal  matrix  with  nonpositive  elements. 


7  Stability  of  the  algorithm 


The  norm  of  an  operator  Ti  is  defined  as 

\\n\\  =  sup  mi  ■  (7.1) 

m=i 

If  the  operator  is  normal,  that  is,  it  commutes  with  its  adjoint,  then  its  norm  coincides 
with  its  spectral  radius  p{Ti),  the  supremum  absolute  value  of  its  eigenvalues.  In  general. 
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p('W)  <  1|H||.  A  family  of  amplification  operators  (matrices)  0At(o!)  is  called  conditionally 
stable  if  there  exists  a  constant  C(t,T)  such  that  [16] 

lie2.(«)ll<c(T,r),  (7.2) 

for  all  At  €  (0,r),  all  0  <  nAt  <  T  for  some  positive  r  and  T,  and  all  parameters  a. 
The  unconditional  stability  of  ^At(«)  means  that  (7.2)  holds  uniformly  in  n  >  0  for  any 
At  >  0  and  for  all  a,  that  is,  C  is  independent  of  T  and  r.  Parameters  a  can  be  all  wave 
vectors  k  supported  by  the  grid  Or  simply  grid  values  of  the  position  vector  x.  They  can 
also  include  parameters  of  the  medium.  Note  that  if  is  not  normal,  then  p(^At)  <  II^Atll 
aad,  therefore,  the  von  Neumann  condition  p{GAt)  <  1  is  no  longer  sufficient  for  stability, 
while  still  being  necessary. 

Theorem  7.1.  For  multiresonant  Lorentz  models,  the  split  algorithm  in  the  field  repre¬ 
sentation  is  unconditionally  stable. 

Proof.  We  shall  prove  that 

WSLW  <  1 ,  (7.3) 

which  leads  to  the  theorem  statement 

ll(6£)”ll  <  lie£ir  <  1  ■  (7-4) 

By  definition  and  making  use  of  the  inequality,  \\AB\\  <  ||A||  ||B|1,  we  get 

\\QL\\  = 

<  (7-5) 

because  is  a  unitary  operator,  so  its  norm  equals  1.  The  operator  e*^  is  block-diagonal 
(see  (4.8)  and  (5.14)).  The  norm  of  a  block-diagonal  operator  is  the  maximal  norm  of  its 
blocks.  The  upper  left  corner  block  is  given  by  the  unitary  operator  whose  norm  equals 
1.  We  have  then 

||e*’^^||  =  max|l,  (7.6) 

The  norm  of  the  exponential  of  cmi  be  found  by  direct  calculation  using  the  fact 
that  ||AlP  =  ||A*A.||  =  p{A*A)  and  the  explicit  form  of  6^^“  given  in  (5.15).  For  small 
attenuation,  uj^  -  yl  -  >  0,  we  define  Zo  =  (ya/i'a)  sm{uat)  so  that  the  largest  eigenvalue 

has  the  form 

||e*WMa|p  = 

=  e-2^“‘(l-b2z^-f-2z„x/r+^)  .  (7.7) 

Since  Za  <  ygt  ^  t/  for  t  >  0,  the  function  (7.7)  is  bounded  from  above  by  f{y)  =  e“^*'(l  -I- 
2y^  -1-  2y  Y^l  +  jp).  It  is  not  hard  to  verify  that  the  derivative  f'{y)  is  negative  for  all  y  >  0, 
and  that  /(O)  =  1.  Hence,  replacing  t  by  At/2,  we  conclude  that 
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from  which  (7.3)  follows.  For  large  attenuation  (like  in  Drude  metals),  —  7a  =  ~K  ^  0, 
in  (7.7)  we  get  Za  =  {la/^a)  sinh(i/at)  =  Za{t).  For  t  >  0  the  latter  relation  defines  the  inverse 
function  t  =  t{za)-  Once  again,  the  derivative  of  (7.7)  with  respect  to  Za  can  be  shown  to  be 
negative  for  all  positive  Za  while  a.t  Za  —  0  the  function  equals  1.  So  inequalities  (7.8)  and 
(7.3)  hold  in  this  case  too.  This  completes  the  proof. 

The  proof  of  Theorem  7.1  given  above  is  not  the  most  economical.  However,  the  idea 
of  estimating  the  norm  of  the  exponential  of  the  matter  Hamiltonian  in  order  to  investigate 
stability  of  the  algorithm  can  be  applied  numerically  to  systems  more  general  than  the 
Lorentz  model  because  Tifj  is  local  on  the  grid,  that  is,  it  does  not  contain  derivatives.  So 
the  exponentials  of  and  its  adjoint  are  not  expensive  to  calculate  numerically  for  some 
trial  values  of  At  to  see  if  (7.8)  holds. 

We  give  an  alternative  proof  of  the  unconditional  stability  in  the  case  of  the  induction 
representation  of  the  multi-resonant  Lorentz  model  where  an  analytical  expression  of  the 
exponent  of  the  matter  Hamiltonian  is  too  hard  to  find,  not  to  mention  its  norm.  We  shall 
make  use  of  the  following  obvious  lenuna. 

Lemma  7.2.  Let  a  vector  ^t,  t  >  0,  be  a  solution  of  the  linear  equation  dti^t  =  (Tf-I- V)^t 
where  the  linear  operators  H  and  V  satisfy  the  conditions  7i*  =  —7i  and  V*  =  V  <  0 
(negative  semidefinite).  Then  l^t|  <  IV’ol  for  all  t  >  0. 

The  proof  follows  from  an  obvious  relation 


(7.9) 

As  a  consequence  we  also  get 

||e*(H+v)||  <  1  ^ 

(7.10) 

for  all  t  >  0. 

Theorem  7.3.  For  multiresonant  Lorentz  models,  the  split  algorithm  in  the  induction 
representation  is  unconditionally  stable. 

Proof.  If  the  attenuation  is  absent,  the  amplification  matrix  ek,  is  unitary  with  respect 
to  the  energy  scalar  product  (/f-scalar  product)  as  is  shown  in  (6.10).  Hence,  ||(^At)”lU  =  1- 
When  the  attenuation  is  switched  on,  the  unitarity  of  the  amplification  matrix  can  get 
violated  only  through  cxp(AtV^)  because  the  operator  cxp(At7f^)  remains  unitary  with 
respect  to  the  //-scalar  product.  The  idea  is  to  prove  the  unconditional  stability  with  respect 
to  the  //-norm.  The  theorem  statement  would  follow  from  the  equivalence  of  the  Euclidean 
and  //-  norms.  Recall  that  two  norms  j'i'l  and  are  equivalent  if  there  exist  two  positive 
constants  Ci  2  such  that 

<  W  <  (721^1  ,  (7.11) 

for  all  All  topological  properties  of  the  space  spanned  by  ^  are  the  same  in  one  norm  as 
in  the  other;  in  particular  convergence  of  a  sequence,  boundedness  of  a  set,  boimdedness  of  a 
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linear  operator,  and  uniform  boundedness  of  a  family  of  operators  are  all  invariant  concepts 
under  a  change  of  one  norm  to  the  other.  If  |^|'  =  then  ||.A|1^  =  and 

m,  =  <  115-^11 1^1 , 

1^1  =  <  IW =  ll^ll  > 

so  that  the  two  norms  are  indeed  equivalent 

||S|r‘|«'l  <  l«l„  <  l|5“‘||  l^l  .  (7.12) 

Since  ||>1W||  =  ||>l|l  for  any  bounded  operator  A  and  a  unitary  operator  U,  we  infer  that 

ii(6i.riu  <  iieiiic = ■  (7.13) 

When  7a  =  0  (no  attenuation),  the  operator  is  skew-symmetric  (cf.  (6.9)).  When 

7„  ^  0,  the  operator  S~^V^S  acquires  an  addition  which  is  a  diagonal  operator  with  nonpos- 
itive  elements  as  follows  from  (2.25)  and  (2.28).  Therefore  the  inequality  (7.10)  must  hold 
for  it  as  a  consequence  of  Lemma  7.2, 

from  which  the  uniform  boundedness  of  the  family  (Git)"  respect  to  the  //-norm  im¬ 
mediately  follows.  By  the  equivalence  of  the  two  norms  (7.12),  the  family  (Git)”  is  also 
uniformly  boimded  in  the  Euclidean  norm, 

ll(ei,)”ll<l|S|ll|S-'l|,  (7.14) 

which  completes  the  proof. 

Comment.  The  same  idea  of  making  use  of  the  norm  equivalence,  which  actually  goes 
in  line  with  the  Kreiss  matrix  theorem  (its  last  part)  [23,  16],  can  be  applied  to  analyze  the 
stability  of  the  split  algorithm  for  generic  passive  media.  It  is  not  hard  to  find  a  quadratic 
T . fl.gr angi  an  local  in  time  such  that  the  corresponding  Euler  equations  describe  propagation 
of  an  electromagnetic  pulse  in  generic  non-absorbing  media.  Due  to  time  translation  symme¬ 
try,  the  system  should  have  a  conserved  quantity  according  to  the  Noether  theorem  [8].  This 
integral  of  motion  coincides  with  the  canonical  Hamiltonian  which  is  a  quadratic  form  of 
if  the  linear  response  approximation  is  valid.  By  analogy  with  the  //-norm,  one  could  try 
to  identify  the  canonical  Hamiltonian  with  the  new  norm  of  which  is  conserved  by  con¬ 
struction  and,  hence,  in  an  attenuation-free  medium  the  corresponding  evolution  operator 
is  unitary  Thus,  it  would  always  be  possible  to  arrange  the  split  so  that  the  amplifica¬ 
tion  operator  is  unitary  too.  From  the  physical  point  of  view,  it  is  then  naturally  expected 
that,  when  absorption  is  added  to  the  system,  the  attenuation  operator  would  generally 
satisfy  the  condition  (6.11)  because  Fourier  amplitudes  of  fields  are  exponentially  attenu¬ 
ated  in  passive  media.  The  latter  would  make  it  possible  to  apply  Lemma  7.2  to  prove  the 
unconditional  stability  of  the  amplification  operator  with  respect  to  the  norm  defined  by 
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the  canonical  Hamiltonian  along  the  lines  similar  to  the  proof  of  Theorem  7.3.  An  obstar 
cle  for  this  rather  natural  idea  to  generalize  Theorem  7.3  to  generic  passive  media  is  that 
the  canonical  Hamiltonian  is  not,  in  general,  positive  definite.  It  becomes  positive  only  on 
solutions  of  the  equations  of  motion  for  matter  fields,  which  is  a  rather  common  feature  of 
T.fl.grfl.Tigifl.n  systems  with  higher  order  time  derivatives.  Thus,  the  canonical  Hamiltonian 
does  not  always  define  a  positive  definite  quadratic  form  in  the  Hilbert  space  for  a  generic 
passive  media,  and,  hence,  cannot  serve  as  a  new  (conservative)  norm.  The  study  of  con¬ 
ditions  on  attenuation-free  media  imder  which  a  positive  definite  and  conserved  quadratic 
form  does  exist  goes  beyond  the  scope  of  this  paper  since  it  would  require  the  canonical 
formalism  and  the  Noether  theorem  for  theories  with  higher-order  time  derivatives,  which 
is  rather  involved  for  generic  peissive  media.  The  question  can  be  addressed  more  easily  for 
each  particular  medium  model  of  interest.  However,  the  unconditional  stability  might  be 
excessive  as  far  as  practical  needs  are  concerned.  It  is  more  important  to  make  the  split 
algorithm  convergent  for  a  generic  passive  medium.  Then  one  should  use  the  equivalence  of 
(conditional)  stability  and  convergence  according  to  the  fundamental  convergence  theorem 
due  to  Kantorovich  [24,  16]. 

Our  findings  in  this  latter  approach  are  summarized  in  the  following  theorem. 

Theorem  7.4.  Suppose  that  the  medium  response  function  satisfies  the  causality  con¬ 
ditions  (that  is,  its  Fourier  transform  has  poles  only  in  the  lower  half  of  the  frequency  plane, 
hnw  <  0).  Let  lit  be  an  exact  evolution  matrix  in  the  grid  representation  (as  defined  in 
Section  3),  and  ^At  be  an  amplification  matrix  in  some  third  order  split  algorithm.  Then 
for  band- width  limited  wave  packets  the  split  algorithm  is  (conditionally)  stable  and  for  all 
0  <  n  <  N,  T  =  NAt,  and  0  <  At  <  r,  there  exist  a  constant  Cm,  which  depends  only  on 
the  medium  parameters,  and  a  constant  Wm,  which  depends  also  on  r,  such  that 
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<  C„  +  i(M,T), 

(7.15) 

8{At,T) 

=  C„  - 1)  =  O(At^)  . 

(7.16) 

1 

< 

<  5{At,T). 

(7.17) 

Remark.  Before  proving  the  theorem,  let  us  discuss  its  significance  for  practical  ap¬ 
plications.  Inequality  (7.15)  implies  conditional  stability,  while  (7.17)  establishes  a  relation 
between  the  accuracy  (and  convergence)  of  the  split  approximation  and  the  imiform  botmd 
in  the  stabihty  condition  (7.15).  By  making  the  time  step  At  smaller,  any  desired  accuracy 
can  be  achieved  during  the  total  (fixed)  simulation  time  T.  The  latter  impfies,  of  course, 
that  the  grid  is  assumed  to  be  chosen  fine  enough  (in  accord  with  the  Shannon  sampling 
theorem)  to  accurately  reproduce  the  initial  pulse  configuration  via  the  fast  Fourier  method. 
Indeed,  let  =  GAt^o  be  a  simulated  solution,  and  =  Ut'^o  be  an  exact  solution,  then 
from  (7.17)  it  follows  that 

\%  -  <  S{At,T)\%\  =  0(At")  ,  (7.18) 

for  all  0  <  t  <  T  and  any  fixed  total  simulation  time  T  which  is  roughly  ‘iLjc  where  L  is  the 
simulation  box  aze  and  c  the  speed  of  light.  Now  we  turn  to  the  proof. 


23 


Proof.  In  the  grid  Fourier  basis,  =  12k  where  k  spans  the  dual  lattice. 

By  construction  of  the  Hamiltonian,  each  Fourier  mode  (k)  evolves  exactly  as  in  the  con¬ 
tinuum  case.  Since  the  medium  response  function  satisfies  the  causality  conditions,  Fourier 
amplitudes  of  the  electromagnetic  and  response  fields  as  well  as  their  time  derivatives  are 
bounded  functions  of  time.  The  amplitudes  cannot  grow  infinitely  large  because  of  dissipa¬ 
tion  [19].  The  number  of  Fourier  modes  is  finite  on  the  grid  (only  bandwidth  limited  initial 
wave  packets  are  considered)  and,  hence,  j'l'f  |  <  Cq  for  all  t  >  0  because  components  of  the 
auxiliary  field  8^6  linear  combinations  of  the  response  field  and  its  time  derivatives.  The 
latter  inequality  is  equivalent  to  the  evolution  matrix  being  uniformly  bounded  for  alH  >  0, 

\\llt\\<Cm.  (7.19) 

Let  UAt  -  OAt  =  and  Wa*  =  Wo  +  0(At)  for  small  At  according  to  a  third  order 

split  (cf.  (4.2)  -  (4.4)).  Let  Wm  ==  CUsup^t  ||WAt|l  for  0  <  At  <  r  and  some  positive  finite 
r.  Using  the  semigroup  property  =  Z4a<  and  (7.19)  we  infer  that 


11^.11  - 

\K,  -  M.  -  6S.)II 

< 

IM  +  IK.-ail 

(7.20) 

< 

1  \ 

(7.21) 

=r 

/  71—X 

Cm+\\-AtUY^UAt(n-k- 

V  fc— n 

-l)WAtZ4At  j  +  •  •  •  II 

(7.22) 

< 

\  w — u 

Cm  +  Cm[{l  +  At^Wm)^- 

1)] 

(7.23) 

< 

Cm  +  S{At,T). 

(7.24) 

Inequality  (7.17)  readily  follows  from  comparing  the  right  hand  side  of  (7.20)  with  those  of 
(7.21)-(7.24).  This  completes  the  proof. 


8  Convergence  and  accuracy  analysis 


To  estimate  the  accuracy  of  the  algorithm  at  a  fixed  finite  grid  size  N,  consider  the  following 
quantity 

UN,  Ai)  =  I  -  {5i,r)  f?!  /  i®oi  <  ||(«f,)“  -  (e2)”||  (S.1) 


which  specifies  a  deviation  of  the  approximate  solution  from  the  exact  one  relative  to  a 
given  norm.  Here  U^t  is  an  exact  evolution  operator.  The  accuracy  estimate  l3n{N,  At)  is  a 
norm  dependent  quantity.  The  choice  of  norm  is  usually  determined  by  practical  needs.  We 
use  the  norm  related  to  the  electromagnetic  energy  of  the  system  and  investigate,  first,  the 
behavior  of  /3n(^)  At)  as  At  goes  to  zero,  while  Atn  —  t  remains  fixed  and  does  not  exceed 
some  positive  constant,  t<T. 
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Theorem  8.1.  For  multi-resonaat  Lorentz  models,  there  exists  a  positive  constant 
such  that 

<At^TW%,  (8.2) 

where  cf  =  1  and  cj  =  |15|1 1|5“^||,  for  all  At  e  (0,  r)  and  nAt  <T. 

Proof.  In  the  field  representation  Q  —  F,  ||(W|^)"||  <  1  and  ||^At)"||  <  1  for  any 
integer  n,  as  a  consequence  of  Lemma  7.2  for  the  multi-resonant  Lorentz  model.  The  same 
inequalities  hold  in  the  induction  representation  if  the  norm  is  replaced  by  the  ^-norm. 
According  to  the  split  algorithm  (4.2)-(4.4),  =  At^W^i-  Let  =  sup^t  llVVAtH 

for  At  €  (0,  r)  for  some  positive  r  (a  maximal  time  step  used  in  simulations).  We  then  have 
the  following  chain  of  inequalities  that  lead  to  the  theorem  statement 

||(w£)”-(®"||  =  (wl.  -  ea  (e£)”“‘  (8.3) 

A:=l 

<  (n-l)At^l|>Vl,||  (8.4) 

<  At'^TW^  ,  (8.5) 

In  the  case  of  the  induction  representation,  inequality  (8.4)  holds  relative  to  the  ^-norm.  The 
theorem  statement  (8.2)  follows  from  the  norm  equivalence  (7.12),  c7^|jA||  <  ||.4||^  <  c/|1.4|| 
for  any  operator  A.  The  proof  is  complete. 

Remark.  In  simulations,  the  continuum  limit  AT  — oo  is  never  achieved.  Hence  the 
operators  in  the  split  algorithm  (4.1)  remain  bounded  versus  the  unboimded  case  of  (4.1). 
It  is  known  that  the  convergence  rate  of  Ai(oo,  At)  as  At  0  estimated  by  the  operator 
norm  as  in  the  right  hand  side  of  (8.1)  is  no  longer  of  order  O(At^)  but  rather  of  0{y/M) 
(see,  e.g.,  [25]  and  references  therein).  For  unboimded  operators,  the  estimate  (8.5)  is  not 
valid.  This  suggests  that  the  convergence  rate  /3n(A,  At)  may  depend,  even  sigmfic^tly,  on 
the  initial  vector  'io  bs  N  increases. 

In  a  general  case,  the  quantity  5(At,T)  in  Theorem  7.4  determines  the  accuracy  of  the 
split  algorithm  with  respect  to  the  norm  (7.1)  on  a  finite  grid.  To  make  simulation  errors 
small,  it  is  sufficient  to  require  that 

l|WAt-a2tl|  =  Ai"|l>V^,||«l.  (8.6) 

Making  use  of  (4.4)  and  the  fact  that  the  norm  of  a  matrix  does  not  exceed  the  maximal 
norm  of  its  blocks,  we  infer  for  a  multi-resonant  Lorentz  model  that,  in  order  for  (8.6)  to 
hold,  the  following  inequalities  are  sufficient: 

CJpaAt  <5C  1  ,  (jJmaxAt  «:  1  ,  WaAt  <K  1  ,  7aAt  «  1  ,  (8.7) 

and,  yet  another  one, 

cAt  «  1  .  (8.8) 
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Here  Wmax  is  the  Tnayimal  frequency  of  the  initial  wave  packet.  The  right  hand  side  of  (8.6) 
is  a  sum  of  two  types  of  terms.  There  are  terms  of  the  cubic  order  in  numbers  (8.7)  as  well 
as  a  term  Unear  in  (8.8)  with  the  coefficient  being  quadratic  in  (8.7).  The  ratio  in  (8.8)  can 
roughly  be  estimated  from  IVwpa]  <  a>pa/Ar  with  Ar  being  the  grid  step.  The  condition 
(8.8)  impUes  then  that  the  distance  traveled  by  the  wave  packet  during  one  time  step  should 
be  much  smaller  than  the  grid  step. 

To  complete  the  discussion,  one  should  also  analyze  the  accuracy  of  the  Gauss  law  (2.3). 
Note  that  the  constraints  are  automatically  fulfilled  in  the  continuum  theory  due  to  the 
Dirac  involution  relations  (1.3).  By  projecting  the  continuum  theory  onto  a  finite  grid 
and  replacing  the  exact  evolution  operator  by  its  approximation  in  the  split  algorithm,  the 
involution  relations  might  be  violated,  thus  leading  to  errors  and  potential  instabifities  of 
the  algorithm.  A  good  example  of  this  kind  is  numerical  general  relativity  (although  the 
nonlinearity  of  the  equations  of  motion  plays  the  central  role  in  generating  instabilities  due 
to  the  violation  of  the  Dirac  involution  relations). 

It  is  not  hard  to  be  convinced  that  the  Gauss  law  (2.3)  is  equivalent  to  the  following 
constraint  on  state  vectors 

C'»,'  =  0,  C'=(^  5),  C=(J  J)?^||,  =  (8.9) 

The  operator  Vw  projects  a  vector  field  onto  its  longitudinal  component.  In  other  words, 
it  acts  as  the  identity  operator  if  the  vector  field  is  conservative,  and  it  annihilates  any 
rotational  vector  field  (which  is  the  cmrl  of  another  vector  field).  In  the  field  representation 
we  get  =  S~^C^S  with  S  defined  in  (2.21).  On  the  grid,  the  action  of  the  operator 

is  defined  by  the  rule  (3.5),  that  is,  by  (8.9)  we  understand  —  0.  Thus, 

the  Gauss  law  requires  that  the  Fourier  transform  of  the  inductions  should  not  acquire 
components  parallel  to  wave  vectors  of  the  dual  grid.  This  is  obviously  guaranteed  if  the 
exact  evolution  operator,  Vt^  —  exp(t7i^),  is  used  to  generate  the  time  evolution  because 

qQ^Q  =  0  ,  (8.10) 

and,  hence,  =  0.  A  problem  may  arise  when  the  approximate  evolution 

operator,  (^^j)”,  is  used  to  evolve  the  initial  wave  packet  Prom  linearity  of  the  system, 
it  is  natural  to  expect  that  the  Gauss  law  violation  should  be  of  the  same  order  as  the 
accuracy  of  a  raunerical  solution  of  dynamical  Maxwell’s  equations.  However,  we  shall  take 
a  closer  look  at  the  problem  and  find  a  pleasant  result  important  in  practice,  which  is  stated 
in  the  following  theorem. 

Theorem  8.2.  Assuming  linear  response  theory  for  any  passive  medium,  the  Gauss  law 
holds  exactly  in  the  split  algorithm  in  the  induction  representation. 

Proof.  In  the  induction  representation,  identity  (8.10)  is  equivalent  to  two  identities  for 
the  blocks  of  C^VP,  namely,  CHo  =  0  and  CVjm  —  0.  The  first  one  is  obvious.  The  second 
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one  follows  from  (2.24)  established  for  any  passive  medium.  The  key  observation  is  that  the 
identity 

=  0  (8.11) 

holds  thanks  to  the  two  above  identities  and  (4.11).  Indeed,  in  the  Fomier  basis  (8.11)  is 
equivalent  to  the  vanishing  of  the  triple  vector  product  k  •  (k  x  A)  for  some  A  regular  at 
k  =  0.  Then  from  (8.9)  and  (8.11)  it  Mows  that 

-7i^)  =  Q.  (8.12) 

As  a  consequence  of  (8.11)  and  (8.12),  we  infer  that 

{GitT  H  ^  ^  ^  Q  ^  ^3^ 

which  is  the  statement  of  the  theorem. 

In  the  field  representation  the  Gauss  law  can  be  enforced  by  means  of  the  projection 
formalism  discussed  in  Section  1.  The  projection  operator  is,  obviously,  V  =  1-C^.  Its 
action  is  computed  in  the  grid  representation  by  the  fast  Fourier  method  according  to  the 
rule  (3.5).  Without  the  use  of  the  projection  formalism,  the  acciuracy  of  the  Gauss  law  is 
stated  in  the  following  technical  proposition. 

Proposition  8.3.  Let  W  =  llWfJI  and  Wc  =  UpMVjtjUC'^l  +  (5(At,T))  where  Cm 
and  5{At,  T)  axe  defined  in  Theorem  7.4,  then 

\C^  {GLT  I  /Kl  <  TWcAt^  +  At^WWce^'^’^*\A^  +  T^/2)  =  0{A^)  ,  (8.14) 

for  all  0  <  n  <  A^,  T  =  NAt  and  any  positive  At. 

Proof.  Since  =  C^,  assuming  that  the  initial  state  satisfies  the  Gauss  law  we 
get 

c^(.sLr^S  =  -c'' {(«£)“ -{e£)"}< 

k-l 

=  .(8.15) 

fc=l  ■  fc=:l 

Denoting  the  left  hand  side  of  (8.14)  by  we  infer  from  (8.15),  by  taking  the  norm  of  both 
sides,  that  ^ 

a„  <  IIIC",  Wiilll  X!  ll(eA.r‘ll  +  A«^|K,||  , 

fc=l  k=l 

forn  >  1  and  ai  <  WcAf.  By  Theorem  7.4,  powers  of  the  amplification  matrix  Gm  are 
bounded.  Hence  the  following  recursion  inequality  holds 

a„  <  (n  -  l)At^Wc  +  At^W{an-i  +  a„_2  +  •  •  •  +  Oi)  .  (8.16) 
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Iterating  (8.16)  n  —  1  times,  we  deduce  that 

On  <  {n-l)Ai^Wc  +  AfW{{n-2)AfWc  +  il  +  ^1^W){an-2  +  an-3  +  ---ai)} 

{n-2 

Y,in  -  2  -  A;)(l  +  Al^Wf  +  (1  +  At^Wf-^ 

k=0 

One  can  find  an  explicit  form  for  the  sum  in  the  latter  equation.  However,  it  is  a  cumbersome 
expression.  For  practical  purposes,  we  give  a  simpler  estimate  which  is  stated  in  (8.14).  First, 
factor  out  (l+At®VF)”“^  in  the  brackets,  and  then  use  obvious  inequalities  (l+At®VF)~*'  <  1 
and  (1  +  <  exp{TWA^),  which  leads  to  (8.14). 

In  the  case  of  the  Lorentz  model,  Wc  =  ||  [C^,  Wit]  II  because  all  powers  of  the  amphfica- 
tion  matrix  are  uniformly  bounded  by  1.  For  small  At,  a  good  estimate  can  be  obtained  by 
computing  Wc  for  At  =  0  using  (4.4). 

The  convergence  rate  as  the  number  of  grid  points  N  increases  is  determined  by  the  con¬ 
vergence  rate  of  the  fast  Fourier  transform  which  is  exponential  versus  polynomial  in  finite 
difference  schemes,  provided  parameters  of  the  medium  are  smooth  functions  of  position 
[1,  5].  As  is  well  known  fi-om  Fourier  analysis,  the  convergence  rate  can  be  affected  for  func¬ 
tions  which  have  discontinuities  [5].  The  latter  is,  unfortunately,  the  case  in  electromagnetic 
scattering  problems.  Suppose  there  is  an  interface  between  two  media.  It  can  be  deduced 
fi:om  the  Maxwell’s  equations  that  the  components  of  the  electric  and  magnetic  fields,  Et 
and  Ht,  tangential  to  the  interface  miist  be  continuous,  provided  there  is  no  surface  electric 
current  on  the  interface.  Prom  the  Gauss  law  it  follows  that  the  components  of  the  induc¬ 
tions,  Dt  and  Bt,  normal  to  the  interface  must  be  continuous,  provided  there  is  no  surface 
charge  on  the  interface.  In  contrast,  the  normal  components  of  the  fields  and  the  tangential 
components  of  the  inductions  can  be  discontinuous.  Their  discontinuities  are  proportional  to 
discontinuities  of  medium  parameters  (e.g.,  discontinuities  in  plasma  frequencies  in  Lorentz 
models).  Therefore,  in  either  the  induction  or  field  representation,  there  are  components 
which  suffer  discontinuities  at  the  interface.  Consequently,  the  convergence  rate  of  the  split 
algorithm  for  Maxwell’s  theory  might  be  slower  than  that  in  quantum  mechanics  with  a 
discontinuous  potential  because  in  the  latter  case  the  wave  function  remains  continuous. 

Another  source  of  errors  that  affects  the  convergence  rate  as  N  increases  is  the  aliasing 
problem  in  the  fast  Fourier  transform  on  the  grid.  Note  that,  even  though  the  initial  wave 
packet  is  band-width  limited  and  the  grid  is  chosen  fine  enough  to  eliminate  errors  in  doing 
its  fast  Fourier  transform  back  and  forth,  the  wave  packet  looses  this  property  after  the 
operator  e3qp(AfV®)  is  apphed  to  it.  As  a  result,  the  aliasing  problem  arises  in  spatial 
domains  where  varies  (typically  at  interfaces  between  different  types  of  media). 

The  above  two  problems  that  also  reduce  the  accuracy  of  the  algorithm  are  well  known 
and  studied  in  the  theory  of  the  fast  Fourier  transform  [5].  The  only  way  to  cope  with 
them  is  to  make  the  grid  finer  in  the  areas  where  medium  parameters  have  discontinuities. 
However,  the  fast  Fourier  algorithm  requires  a  uniform  equispaced  lattice,  which  might  lead 
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to  wasting  computer  resources  if  the  increased  resolution  is  necessary  only  in  relatively  small 
areas  of  the  computational  volume  of  the  problem  (e.g.,  only  near  an  interface  between  two 
media).  There  are  several  ways  to  modify  the  algorithm  when  the  above  problems  are  too 
expensive  to  overcome  by  making  a  uniform  grid  finer. 

First,  the  grid  can  be  made  fine  enough  so  that  the  action  of  powers  of  the  Hamiltonian 
'HP  on  the  state  vector  is  sufficiently  accurate  in  the  Fourier  basis  as  specified  by  the 
rul^  (3.4)  and  (3.5).  The  operator  is  projected  onto  the  Krylov  space  spanned  by 
vectors  (7f‘2)*^^o>  *  =  0,l,...,n.  Its  exponent  (the  evolution  operator)  is  then  computed 
by  diagonalizing  "H^  instead  of  using  the  Lie-Trotter  formula.  Usually,  it  is  sufficient  to 
take  a  low  dimensional  Krylov  space.  This  method  is  known  as  the  Lanczos  method  [26].  A 
detailed  study  of  the  Krylov-Lanczos  method  as  well  as  other  similar  pseudospectral  methods 
in  Maxwell  theory  will  be  given  elsewhere. 

Second,  one  can  give  up  a  uniform  grid,  while  preserving  basic  advantages  of  pseudospec¬ 
tral  methods  such  as,  e.g.,  exponential  convergence.  A  possible  way  to  emulate  a  non-umform 
grid  in  a  multiscale  problem  is  to  use  wavelet  bases.  The  problem  here  is  to  compute  the 
action  of  exp(AtW^)  in  the  split  algorithm  because  the  derivative  operator  V  is  not  diagonal 
in  this  basis  (in  contrast  to  the  Fourier  basis).  However,  is  expected  to  be  sparse  in  a 
wavelet  basis  so  that  its  direct  diagonahzation  might  not  be  expensive,  and  a  significant 
reduction  of  computational  costs  can  be  achieved  in  the  spht  algorithm,  by  using  the  fast 
wavelet  transform,  as  compared  to  that  in  the  Fourier  basis.  Otherwise,  the  use  of  (4.10) 
might  be  helpful  in  place  of  the  direct  diagonalization  method.  This  approach  has  proved  to 
be  successful  in  solving  multiscale  initial  value  problems  for  the  Schrodinger  equation  [27]. 
In  the  framework  of  Maxwell’s  theory  for  passive  media,  additional  studies  of  several  issues 
in  time  domain  wavelet  based  algorithms,  like,  e.g.,  stability,  would  still  be  needed. 

Third,  the  fast  Fourier  transform  algorithm  remains  in  place  but  is  applied  to  an  auxiliary 
uniform  grid  that  is  related  to  a  non-uniform  grid  in  physical  coordinates  by  a  change  of 
variables.  Consider  a  change  of  variables  y  =  y(x).  A  uniform  grid  in  the  new  variables 
y  would  generate  a  non-uhiform  grid  in  the  original  Euclidean  (physical)  coordinates  x.  A 
desired  local  density  of  grid  points  in  the  physical  space,  to  enhance  the  sampling  efficiency 
in  designated  regions,  can  be  achieved  by  an  appropriate  choice  of  the  functions  y(x)  [28].  By 
necessity,  the  auxiliary  grid  spans  a  rectangular  box  (with  periodic  boundary  conditions).  Its 
pre-image  in  the  physical  space  would  not  be  a  box  in  general,  save  for  the  case  when  the  map 
y(x)  splits  into  three  individual  one-dimensional  maps  yj  =  Vjixj),  j  =  1,2,3.  Since,  the 
derivatives  are  transformed  as  Vx  =  •A(y)Vy  where  the  3x3  matrix  A  is  position  dependent, 
the  operator  cannot  be  kept  in  the  exponential.  The  action  of  its  exponential  on  the 
state  vector  can  be  approximated  by  the  leapfrog  method  in  which  only  the  action  of  7i^  on 
is  required.  The  latter  can  be  done  by  the  fast  Fourier  method  according  to  the  rules 
(3.5)  and  (3.4)  applied  to  an  operator  being  a  product  of  position  and  derivative  dependent 
operators.  In  contrast  to  the  well  studied  quantmn  mechanical  case,  the  algorithm  appears 
to  be  unstable  for  media  with  absorption.  In  Section  9  a  modification  of  the  leapfrog  scheme 
is  proposed  to  achieve  conditional  stability. 
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The  temporal  leapfrog  scheme 


Here  we  discuss  a  temporal  finite  difference  scheme  applied  to  the  Maxwell  theory  for  passive 
media  in  the  Hamiltonian  formalism.  As  has  been  pointed  out,  such  a  scheme  might  be 
helpful  for  reducing  computational  costs  by  using  non-uniform  grids  in  combination  with 
some  pseudospectral  methods  (e.g.,  wavelet  bases  or  the  fast  Fourier  method  with  a  change 
of  variables).  A  temporal  finite  difference  scheme  can  be  obtained  by  the  following  procedure. 
Let  us  integrate  (1.1)  over  the  interval  {t,t  +  nAt).  We  have 

t+nAt  -V 

%+nAt  =  %  +  n  I  dT^r  =%  +  AtW  %t+fcAt  j  +  0(Ar+l)  ,  (9.1) 

where  the  coefficients  used  to  approximate  the  integral  are  well  known  for  any  n  as 
well  as  the  accuracy  of  the  approximation.  For  example,  one  can  use  the  3/8  Simpson  rule 
for  n  =  3  or  Bode’s  rule  for  n  =  4.  The  iterating  scheme  allows  one  to  compute  the  wave 
function  at  the  sequential  moment  of  time  if  it  is  known  for  n  preceding  moments  of  time. 
Only  the  simplest  case  n  =  2,  for  which  the  mid-point  approximation  for  the  integral  is 
taken,  leading  to  =  0  and  =  2,  will  be  considered  in  detail.  It  is  also  known  as  the 
leapfrog  scheme; 

=  ^t-At  +  •  (^-2) 

The  action  of  the  Hamiltonian  is  computed  in  a  suitable  basis  (as  has  been  noted  above). 
Apart  from  violation  of  the  dispersion  relation  of  electromagnetic  waves,  temporal  finite 
difference  schemes  would  generally  be  unstable  in  media  with  absorption,  in  contrast  to  the 
quantum  mechanical  case.  The  reason  is  that  the  Hamiltonian  in  (9.2)  is  not  anti-Hemutian. 
Consequently,  convergence  to  the  continuum  solution  would  also  be  violated. 


A  general  solution  to  (9.2)  can  be  written  in  the  form 

ai?  =  HAt±  +  , 


(9.3) 

(9.4) 


for  some  initial  state  vectors  'i'o  and  I" At  ( the  vectors  are  determined  by  them).  Stability 
requires  that  there  exists  a  positive  constant  C  such  that 


|^nAt|<C'(|'^+!-f-|'F_|)  ,  (9.5) 

for  all  0  <  n  <  iV,  T  =  NAt  and  0  <  At  <  t.  Note  that  in  general  a  solution  of  (1.1) 
may  have  a  legitimate  exponential  growth  if  the  hermitian  part  of  the  Hamiltonian,  H+H*, 
is  not  negative  semidefinite  (see  Lemma  7.2).  For  this  reason,  a  typical  stability  criterion 
would  be  equivalent  to  the  condition  [16]  that  there  exists  some  positive  constant  Ki  such 
that  11^^^  II  <1-1-  KiAt  uniformly  for  all  parameters  of  G^t  for  0  <  At  <  r,  which  is 
clearly  the  case  for  (9.4)  if  W  is  bounded.  The  latter  leads  to 
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so  that  a  legitimate  exponential  growth  of  the  solution  is  allowed,  i.e.,  C  ~  in  (9.5). 
For  passive  media,  the  Hamiltonian  satisfies  the  conditions  of  Lemma  7.2  and,  hence,  no 
legitimate  exponential  growth  should  be  present  in  a  numerical  solution  in  order  to  achieve 
convergence.  However,  as  we  shall  see  shortly,  the  scheme  (9.2)  always  generates  an  expo¬ 
nentially  growing  solution  for  media  with  attenuation. 

Let  complex  numbers  z  =  be  eigenvalues  of  W  Af.  Since  the  spectral  radius  of 

does  not  exceed  its  norm,  the  necessary  (von  Neumann)  condition  to  suppress  an  exponential 
growth  of  the  solution  (9.3)  reads 


p{gS)  = 


z  ±  vT+l^ 


<  1 


(9.7) 


The  aim  is  to  analyze  the  domain  D  of  the  complex  plane  for  which  (9.7)  holds.  Let 
T]  =  VTT^  and  |7/|  =  r.  The  two  inequalities  in  (9.7)  require  that  for  z  E  D,  <  1, 

where  ^  =  ZT)+zfj.  By  combining  the  latter  inequalities,  one  gets  R^+r^  <  1  or  <  (1— J?^)^. 
On  the  other  hand,  r*  =  1  +  R^  +  2R^  cos{2(p).  Hence,  cos(2(^)  <  -1  which  is  only  possible 
if  ^  =  ±7r/2  or  z  =  ±ii?.  The  necessary  condition  (9.7)  is  satisfied  if 


if  =  ±7r/2  ,  R^  <  1  . 


(9.8) 


This  does  not  yet  guarantee  that  there  is  no  norm  growth.  A  norm  growth,  which  is  poly¬ 
nomial  in  time,  can  still  occur. 


Let  us  investigate  general  properties  of  the  solution  of  (1.1)  when  the  Hamiltonian  satisfies 
the  von  Neumann  condition  (9.8).  For  any  matrix  H  there  exists  a  similarity  transformation 
so  that  S~^HS  has  the  Jordan  normal  form.  Let  be  a  block  of  the  Jordan  normal  form 
corresponding  to  an  eigenvalue  z  of  Tf.  Any  block  hk  is  a,  q^x  bi-diagonal  matrix,  5^  >  1, 
with  all  the  elements  of  the  diagonal  being  equal  to  z  and  all  the  elements  on  the  upper 
superdiagonal  being  equal  to  one.  For  q^  —  z  is  just  a  complex  number.  The  norm 

of  any  solution  of  (1.1)  cannot  grow  faster  than  ||  exp(t7f)||.  Let  a  gz-dimensional  vector 
satisfy  the  equation  dt(i>t  =  Kc/h-  For  a  generic  initial  condition,  the  solution  norm  grows 
polynomially,  |^t|  =  ast  oo.  Using  the  similarity  transformation  5,  we  define 

the  corresponding  //-norm  of  state  vectors  and  the  corresponding  matrix  norm  (cf.  (6.5))  by 
setting  fJL  =  Prom  the  equivalence  of  the  norms  ||  •  ||  and  ||  •  ||p  (see  the  proof  of 

Theorem  7.3),  the  norm  growth  cannot  be  faster  than 

=  max||e*'‘^|l  =  0{P~^)  ,  q  =  raax.q^  ,  t-^oo  ,  (9.9) 

Z  Z 


provided  z  =  ±iR.  However,  a  state  vector  norm  growing  polynomially  with  time  is  unac¬ 
ceptable  from  the  physical  point  of  view  because  in  any  passive  medium  there  is  no  physical 
mechanism  for  such  amplification  of  the  field  amplitudes  in  the  large  time  limit.  Conse¬ 
quently,  we  demand  that  any  model  Hamiltonian  for  a  passive  medium  should  be  similar  to 
a  diagonal  matrix  (i.e.,  7i  is  diagonalizable).  In  this  latter  case,  the  blocks  hz  of  the  Jordan 
normal  form  of  Ji  are  just  complex  numbers  z.  Hence  ||  exp(t/i2)||  =  |exp(±*fi?)|  =  1  so 
that  the  /t-norm  of  any  solution  of  (1.1)  is  conserved  according  to  (9.9). 
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Two  important  conclusions  about  the  leapfrog  scheme  (9.2)  follow  from  our  analysis. 
First,  the  von  Neumann  condition  (9.8)  is  also  sufficient  for  stabiUty.  Indeed,  if  (9.8)  holds 
then  =  p(ai?)  =  1  and,  hence,  l|(ai?)”lU  ^  1  uniformly  in  n  >  0. 

By  the  norm  equivalence,  |1(^^^)"11  is  also  bounded  uniformly  in  n  >  0.  Second,  reversing 
the  argument,  we  conclude  from  the  norm  conservation  of  the  stable  leapfrog  solution  that 
no  attenuation  can  be  added  to  the  Hamiltonian  without  destroying  the  stability  and,  con¬ 
sequently,  the  convergence  to  the  continniun  solution.  Whenever  the  attenuation  is  added, 
the  leapfrog  solution  would  always  contain  an  exponentijdly  growing  component,  while  this 
would  not  be  so  for  a  continuum  solution  (see  Lemma  7.2). 

Since  =  1,  only  one  of  the  two  independent  solutions  in  (9.3)  would  grow  expo¬ 
nentially  whenever  the  attenuation  is  added.  Theoretically,  for  <  0  the  exponentially 

growing  solution  can  be  eliminated  by  choosing  the  initial  condition  so  that  =  0  which 
is  equivalent  to  the  initial  condition  Practically,  this  is  never  possible  due 

to  rounding  errors  and/or  numerical  errors  in  computing  Even  for  a  small  l^'-j  in 

(9.3)  the  growing  part  would  eventually  become  comparable  with  the  exponentially  atten¬ 
uating  solution  generated  by  ^+.  A  r«iuction  of  the  time  step  would  not  be  helpful  since 
the  constant  Ki  in  (9.6)  is  independent  of  At  while  the  simulation  time  T  is  fixed  by  the 
dimension  of  the  simulation  volume  and  the  speed  of  light.  One  needs  at  least  to  modify  the 
scheme  so  that  there  exists  a  constant  Kp  such  that 

||Cli?||<l  +  ApAtP,  p>l,  (9.10) 

for  0  <  At  <  r.  Indeed,  it  follows  from  (9.10)  that  ||(5S^)"I1  <  exp(ifpTAt^“^)  =  1-1- 
0(At*’“^)  for  all  0  <  n  <  W  where  iVAt  =  T.  The  norm  growth  could  be  reduced  as  much 
as  desired  by  making  the  time  step  smaller.  Next  we  show  how  to  modify  the  leapfrog  scheme 
to  make  (9.10)  valid  for  at  least  p  =  3  and,  if  the  Hamiltonian  is  normal,  an  even  stronger 
result  holds,  namely,  ATp  =  0. 

Let  a.  =  'Hq-\-V  where  V*  -f  V  <  0  (negative  semidefinite)  and  Hq  —  —Wo-  In  (1.1)  we 
make  a  substitution  =  exp(tV)^'t.  The  new  state  vector  satisfies  an  equation  with  a 
time  dependent  Hamiltonian, 

dt^t  =  =  Wt  ,  (9-11) 

and  with  the  same  initial  condition  $o  =  ’®'o-  Applying  the  leapfrog  method  to  (9.11)  we  get 

valid  up  to  0(At®).  Returning  to  the  initial  variables,  we  arrive 
at  the  following  recurrence  relation 


’®'t+At  =  iC2At'®'t-At  +  2AtCAtH.o^t  >  (9-12) 

where  =  exp(AtV).  All  the  derivative  operators  are  included  into  the  anti-Hermitian 
part  Wo  of  the  Hamiltonian  H,  while  the  attenuation  operator  V  might  even  be  independent 
of  position  and,  hence,  C^t  has  to  be  computed  only  once  for  given  medium  parameters  and 
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time  step.  It  can  often  be  done  analytically  as,  for  example,  in  multiresonant  Lorentz  models 
(see  Section  10).  On  the  other  hand,  by  Lemma  7.2,  ||£Ai||  <  1  for  any  At  >  0,  and  one 
might  hope  to  stabilize  the  leapfrog  scheme  by  satisfying  the  stability  condition  (9.8)  for  Tio 
only,  that  is,  1  +  TilAt^  is  positive  semidefinite.  This  is  indeed  the  case.  The  amplification 
matrix,  for  the  recurrence  (9.12),  satisfies  the  equation 

GAt  =  C,iAtG^  +  lAtCAtHQ  .  (9.13) 

According  to  our  analysis  of  the  von  Neumann  stability  condition  (9.8),  the  anti-Hetmiticity 
condition  of  Tio  in  (9-12)  and  (9.13)  can  be  weakened  by  demanding  that  Hq  is  related  to 
an  anti-Hermitian  matrix  by  a  similarity  transformation.  Some  important  properties  of  the 
amplification  matrix  obtained  from  (9.13)  are  stated  in  the  following  theorem. 

Theorem  9.1.  Suppose  there  exists  a  similarity  transformation  such  that  S~^TiS  = 
Hs+Vs  where  TiJ  =  —Hs,  the  Hermitian  part  of  Vs  is  negative  semidefinite,  V5+VS  <  0,  and 
Tis  also  satisfies  the  von  Neumann  stability  condition  for  the  leapfrog  scheme,  l+T^At^  > 
0  (positive  semidefinite).  Consider  the  amplification  matrix  GAt  of  the  modified  leapfrog 
scheme  (9.13)  with  7U  =  S'HsS~^  and  V  =  5V5«5“^.  Then  there  exists  a  norm  ||  •  ||/, 
equivalent  to  ||  •  ||  such  that  GAt  has  the  following  properties: 

(A)if[T^o,V]  =  0, 

1154.11.  <  1  PM) 

uniformly  in  n  >  0; 

{B)  if  [Tio,  V]  0)  there  exists  a  non-negative  constant  such  that 

\\GAt\\^<l  +  KiAe  ,  (9.15) 

for  0  <  At  <  r  and  some  positive  r. 

Proof.  Part  (A).  If  Tio  and  V  cosnmute,  the  amplification  matrix  GAt  —  satisfies 

(9.13),  provided  Gai  satisfies  the  same  equation  for  V  =  0  (or  Cai  —  1),  which  one  can 
easily  check  by  substituting  the  solution  into  (9.13).  Consider  the  norm  associated  with  the 
similarity  transformation  S  of  the  Hamiltonian,  ||A^i|/i  =  ||«S“^A4S||.  The  norms  ||  •  ||^  and  ||  •  || 
are  equivalent  (see  the  proof  of  Theorem  7.3).  Since  Tio  satisfies  the  von  Neumann  stability 
condition  and  is  anti-Hermitian  relative  to  the  //  scalar  product,  ||^AtlU  =  P{G%d  —  ^ 
(according  to  the  analysis  after  (9.7)).  By  Lemma  7.2,  ||£At||#i  ^  1>  sJ^d  we  infer  that 
ll^'AtlU  =  ll(>CA*^?At)"IU  ^  ^  1  uniformly  for  n  >  0. 

Part  {B).  Solving  (9.13)  by  the  perturbation  theory  in  At,  it  is  not  hard  to  find  that 

GAt  -  Git  =  .  Git  =  CmA^a^  ,  (9.16) 

where  /Ca*  is  regular  in  the  vicinity  of  Ai  =  0  and  vanishes  whenever  Tio  V  commute.  On 
the  grid,  Tio  and  V  are  bounded  operators.  Hence  we  can  find  a  constant  K3  —  sup  At  II^AtlU 
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for  some  open  interval  0  <  At  <  r.  Making  use  of  the  inequality  H^XtlU  ^  ll'CAt/2ll^  <  1, 
we  find 

II^^AtlU  =  wolt  +  At'/CAtlU  <  1  +  At^/f3  .  (9.17) 

which  completes  the  proof. 

The  norm  deviation  of  the  solution  generated  by  the  modified  leapfrog  scheme  (9.12) 
from  the  stable  solution  generated  by  is  of  order  0{M^)  for  the  entire  simulation  time 
T  and,  hence,  by  reducing  At  a  possible  norm  growth  can  be  suppressed  as  much  as  desired. 
Indeed, 

lia  -  (6A.)1t  =  llEes‘(eA.-0(eX<)‘IU 

<  At^ir^  110X7'' lU  ^  =  0(At2) 

Since  in  the  continuum  limit  At  — 0,  both  the  amplification  matrices  Q^^  and  Oa*  generate 
the  same  solution  and  all  the  powers  of  the  former  are  imiformly  bounded  by  construction, 
a  natural  question  to  ask  is  whether  one  can  find  a  recurrence  relation  for  the  function 
^Xa<  =  (OXt)”’®'©  which  could  be  used  in  place  of  (9.12).  It  is  not  difficult  to  derive  an 
equation  for  OXt  similar  to  (9.13),  but,  unfortunately,  this  equation  cannot  be  converted  into 
a  simple  recurrence  relation  for  the  wave  function  itself,  like  (9.12),  suitable  for  numerical 
applications. 

It  should  be  noted  that  if  the  operator  Aa*  in  the  modified  leapfrog  scheme  (9.12)  is 
replaced  by  another  such  that  Z^At-Af  =  P(At^)  and  IjAXtlU  ^  convergence 

is  not  violated  because  Part  B  of  Theorem  9.1  still  holds.  This  observation  is  useful  for 
analytic  computation  of  Z^At-  For  example,  in  the  conditions  of  Theorem  9.1,  put  «S  =  1. 
Let  V  =  Vi  +  V2  so  that  both  Vi,2  have  their  hermitian  parts  negative  semidefinite.  By  using 
the  split  (4.2)  we  get 

=  e^‘v,/2gAmgAm/2  +  0(At^)  =  C%t  +  0{At^)  .  (9.18) 

By  Lemma  7.2,  ||Z^Atll  ^  1  At  >  0.  The  operators  Vi,2  can  be  chosen  so  that  their 
exponentials  can  be  computed  analytically. 


10  Examples  of  the  temporal  leapfrog  algorithm 


There  are  many  possibilities  to  split  the  original  Hamiltonian  H  into  two  parts  that  satisfy  the 
conditions  of  Theorem  9.1  and  thereby  to  make  the  leapfrog  scheme  stable  and  convergent. 
Basic  guide  lines  for  doing  that  are  as  follows.  The  Hamiltonian  TZo  should  contain  all 
the  derivative  operators  in  Ti  and,  yet,  the  von  Neumaim  condition  is  easy  to  estabfish  for 
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Ho-  It  would  also  be  helpful  to  have  an  analytic  expression  for  Cm  at  least  up  to  order 
O(At^).  As  an  illustration,  we  discuss  multiresonant  Lorentz  models  and  geometric  optics. 
To  distinguish  between  the  splits  of  the  Hamiltonian  in  the  spht  and  leapfrog  algorithms, 
we  shall  use  an  index  I  (“leapfrog”)  in  the  latter. 


I  Lorentz  models 


In  the  field  representation  of  the  Hamiltonian  for  multiresonant  Lorentz  models,  we  make 
the  following  decomposition 


/  Ho 
V  Vmf 


(10.1) 


Thanks  to  (2.27)  and  HI  =  -Ho,  the  operator  H£  is  anti-Hermitian.  Prom  (2.28)  it  follows 
that  the  Hermitian  part  of  Vf  is  negative  semidefinite  (7a  >  0).  The  exponential  of  Vf 
is  easily  computed  according  to  (5.14)  and  (5.15).  Let  denote  a  six-component  column 
whose  three  upper  components  coincide  with  and  three  lower  components  equal  ^ 
(see  (2.9)-(2.11)).  As  a  result  we  arrive  at  the  following  scheme 


a 

^r-At  +  2At  Vmfc  - 

Stability  is  ensured  if  Hq  satisfies  the  von  Neumann  condition  (9.8). 
satisfy  the  equation 


(10.2) 

(10.3) 

Eigenvalues  of  Hq 


det(^  —  Hqi)  =  det(z^  —  zHq  —  VfmVmf)  =  0 


(10.4) 


where  the  non-negative  integer  q  depends  on  the  munber  of  resonances  in  the  Lorentz  model. 
Non-zero  eigenvalues  satisfy  the  so-called  pencil  equation  whose  theory  is  well  developed 
and  might  be  useful  for  more  general  models  [29].  Here  we  shall  find  a  simpler  (practical) 
criterion  sufficient  for  (9.8)  to  hold.  Since  the  plasma  frequencies  may  depend  on  position, 
we  apply  the  following  general  idea  [30].  Suppose  we  have  a  finite  difference  scheme  with 
variable  coefficients  in  space.  Consider  a  corresponding  finite  difference  scheme  with  frozen 
coefficients.  It  is  obtained  from  the  original  scheme  by  fixing  the  coefficients  to  particular 
values  everywhere  in  space.  A  finite  difference  scheme  with  variable  coefficients  is  stable 
if  all  the  corresponding  finite  difference  schemes  with  frozen  coefficients  are  stable  [30,  16]. 
So  let  us  fix  the  plasma  frequencies  to  particular  values.  The  spatial  dependence  of  the 
eigenfunctions  for  the  pencil  problem  in  (10.4)  is  given  by  a  harmonic  factor  e3q}(ik  •  x)  and 
the  corresponding  eigenvalues  are  z  =  dri-^/c^^  +  a;|,  where  is  defined  in  (5.16).  Let 
fcmnx  be  the  maximal  norm  of  all  wave  vectors  of  the  initial  wave  packet  and  u;^  be  the 
TTiflYiTnal  value  of  Wp  as  a  function  of  position,  then  a  sufficient  criterion  for  stability  reads 


AtyJ(?k^ax  +  (W^)2  <  1  . 


(10.5) 
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The  scheme  (9.12)  becomes  especially  simple  in  the  case  of  small  attenuation,  7a  <  uJa-  In 
the  complex  representation  of  the  auxihary  fields  (2.30)  (cf.  (2.12))  the  matter  Hamiltonians 
Wm.  are  diagonal  and  the  action  of  its  exponential  is  reduced  to  multiplication  by  a  complex 
munber  (see  Section  5). 

The  stability  condition  (10.5)  can  be  improved  if  one  uses  the  induction  representation 
arranging  the  spht  according  to  (4.11),  that  is,  Hqi  =  'Hq  and  V/  =  V^.  In  this  case  the 
conditions  of  Theorem  9.1  are  met  if  instead  of  (10.5)  we  demand  a  weafer  condition 

^tckmax  <  1  •  (10.6) 

To  prove  this,  we  note  first  that  by  the  similarity  transformation  defined  in  (2.22)  we  get 

Then  (10.6)  is  obviously  the  von  Neumann  stability  condition  for  fig,  while  the  Hermitian 
part  of  V|  is  negative  semidefinite  if  7a  >  0.  The  scheme  is  obtained  from  (9.12)  by  replacing 
H  Hq  and  V  ->  as  defined  in  (4.11).  Since  the  left  hand  side  of  (10.7) 
coincides  with  it  can  also  be  viewed  as  the  leapfrog  scheme  in  the  field  representation 
but  with  the  split  different  from  (10.1).  This  illustrates  the  point  that  the  stability  condition 
of  the  scheme  (9.12)  depends  strongly  on  the  choice  of  Tioi.  The  price  for  a  simpler  stability 
condition  in  the  induction  representation  is  the  lack  of  an  expUcit  form  of  CAt-  However, 
this  problem  can  be  circumvented  by  making  use  of  (9.18).  Indeed,  V|  =  +  Vf  where 

is  defined  in  (4.7).  The  exponentials  of  these  operators  are  computed  in  Section  5.  We 
also  have  ||  exp(AtV^)||  =  1  because  =  -V^  and  ||  exp(AtVf  )||  <  1  by  Lemma  7.2  for 
a  non-negative  At.  We  set 

=  (10.8) 

so  that  ll>Cit||,i  =  <  1.  The  operator  (10.8)  differs  from  =  exp(AtV^) 

=  <S  exp(AfV5)<S“^  by  terms  of  order  0(At®)  and,  hence,  according  to  (9.18),  can  be  used 
in  place  of  CAt  in  the  leapfrog  scheme  without  destroying  its  convergence  and  stability. 


II  Geometric  optics 


Another  simple  example  is  the  case  of  geometric  optics.  For  sake  of  simplicity  we  assume  the 
medium  to  have  no  magnetic  properties.  A  generalization  is  straightforward.  Let  e  =  e(x) 
be  the  dielectric  constant  of  the  medium.  If  the  medium  is  not  isotropic,  then  £  is  symmetric 
positive  definite  3x3  matrix  everywhere  in  space.  We  rewrite  Maxwell’s  equations  in  the 
form 


(10.9) 
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where  the  parentheses  in  (e“^)  mean  that  the  induction  is  first  multiplied  by  e  ^  and  then 
the  curl  of  the  resulting  vector  field  is  computed.  Consider  the  scalar  product 

In  the  grid  representation  of  Section  3,  the  integral  is  replaced  by  the  sum  over  grid  points 
and  Hg  becomes  a  finite  matrix.  The  Hamiltonian  is  anti-Hermitian  with  respect  to  this 
scalar  product,  Hcfi  =  —fiKa-  Therefore  the  corresponding  p  norm  is  preserved  in  the  time 
evolution  generated  by  expitHa),  that  is,  The  electromagnetic  energy 

of  the  wave  packet  is  conserved  because  it  is  proportional  to  the  /x  norm  of  the  initial  state 
vector.  Consequently,  we  expect  that  for  a  sufficiently  small  At  the  original  leapfrog  scheme 
(9.2), 

i4+At  =  +  ^Atnafl ,  (10.11) 

becomes  stable.  To  find  a  liufficient  condition  for  stability,  the  same  idea  of  finite  difference 
schemes  with  frozen  coefficients  can  be  used.  It  obviously  leads  to  a  condition  similar  to 
(10.6), 

Atck^^<l, 

where  k^ax  ^  Ihe  maximal  norm  of  all  lyave  vectors  in  the  medium  which  can  be  estimated 
hy  \/p{£)kmax  with  kmax  being  the  maximal  wave  vector  of  the  initial  pulse  in  vacuum.  The 
spectral  radius  p{e)  is  understood  as  the  maximal  spectral  radius  of  ^(x)  over  x.  If  the 
Fourier  basis  is  used  to  compute  the  derivatives,  the  algorithm  does  not  violate  the  Gauss 
law.  However,  the  algorithm  would  not  conserve  the  p  norm  (or  energy),  rather  a  quantity 
which,  in  many  cases,  approximates  the  energy.  Multiplying  (10.11)  by  V’/  using  the  scalar 
product  (10.10),  we  infer  that 

(’/’f+A*.  i’l )  =  .  V^-At)  =  •  •  •  =  ii’it,  ^i)  .  (10.12) 

By  expanding  the  exponential  in  ~  ^^{At'HG)fl>t  ^  Taylor  series  in  both  sides  of 
(10.12)  and  making  use  of  the  anti-Hermiticity  of  Hg,  we  find  that  the  energy  conservation 
violation  is  of  order  0{At^).  Thus,  it  can  be  made  as  small  as  desired  by  reducing  the  time 
step. 


11  Conclusions 


The  initial  value  problem  in  Maxwell  theory  for  passive  media  has  been  reformulated  in  the 
Hamiltonian  formalism.  The  path  integral  representation  of  the  fundamental  solution  of  the 
Hamiltonian  evolution  equation  has  been  used  to  develop  a  time  domain  numerical  algo¬ 
rithm  for  solving  the  initial  value  problem.  The  algorithm  exhibits  the  main  advantages  of 
pseudospectral  methods  for  solving  differential  equations  such  as  an  exponential  convergence 
(and,  hence,  a  greater  accuracy),  the  absence  of  dispersive  errors  and  numerical  efficiency.  In 
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addition,  the  algorithm  is  unitary,  meaning  that  the  energy  of  the  initial  pulse  is  conserved 
whenever  the  medimn  attenuation  vanishes  (Theorem  6.1).  For  widely  used  multiresonant 
Lorentz  models,  the  algorithm  is  unconditionally  stable  (Theorems  7.1  and  7.3),  and,  for  a 
generic  passive  medium,  conditional  stability  can  always  be  achieved  (Theorem  7.4).  As  the 
time  step  At  goes  to  zero,  the  algorithm  accuracy  is  of  order  O(At^)  (Theorem  8.1).  It  is 
possible  to  increase  the  convergence  rate  (accuracy)  up  to  any  desired  order  0(At”),  n  >  2. 
However,  computational  costs  for  increasing  the  accuracy  in  such  a  way  are  not  necessar¬ 
ily  lower  than  those  for  decreasing  the  time  step  in  the  original  algorithm.  An  important 
advantage  of  the  algorithm  is  that  the  Gauss  law  holds  exactly  in  the  process  of  numerical 
simulations  with  no  extra  computational  cost  (Theorem  8.2). 

A  drawback  of  the  algorithm  is  related  to  well  known  problems  of  the  fast  Fomier  method. 
Namely,  a  slower  rate  of  convergence  for  non-smooth  functions  and  aliasing.  This  might, 
perhaps,  limit  the  advantages  of  the  algorithm  in  some  type  of  scattering  problems  with  com¬ 
plex  target  geometries.  Numerical  tests  are  needed  for  a  quantitative  conclusion.  There  are 
several  pseudospectral  methods  for  approximating  the  fundamental  solution  of  the  Hamil¬ 
tonian  evolution  equation  that  can  help  to  circumvent  this  problem.  We  have  analyzed  one 
of  them  and  formulated  its  stability  criteria  in  the  case  of  general  passive  media  (Theorem 
9.1).  Other  methods  will  be  discussed  elsewhere  as  well  as  the  case  when  radiation  sources 
(antennas)  axe  included. 

It  is  beUeved  that  the  proposed  algorithm  would  be  useful  in  nmnerical  studies  of  elec¬ 
tromagnetic  pulse  propagation  in  passive  media  (e.g.,  foliage,  soil,  etc),  photonic  crystals 
and  devices,  and  also  in  applications  to  scattering  problems  with  targets  made  of  dispersive 
materials. 
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Appendix 


I  Initial  pulse  configurations 

In  principle,  any  field  configuration  can  serve  as  the  initial  configuration.  However  it  is 
often  desired  to  have  an  initial  I)ulse  with  some  specifie  properties  (bandwidth,  polarizar 
tion,  direetion  of  propagation,  etc).  Yet,  the  initial  wave  packet  should  be  built  of  radiation 
(propagating)  electromagnetic  fields.  A  simple  method  based  on  the  fast  Fourier  trans¬ 
form  algorithm  is  given  below  to  obtain  initial  configurations  made  of  radiation  fields  with 
designated  properties. 

A  general  solution  of  the  Maxwell  equations  in  vacuum  can  be  written  in  the  form 

E^r)  =  J  dk  -h  do'\^)  =  Cl+^(r)  -h  Ci’^(r)  (12.1) 

Bt(r)  =  ^  Vx  (cW(r)-Ci->(r))  .  (12.2) 

The  representation  (12.2)  for  the  magnetic  induction  follows  from  the  Maxwell  equations 
and  that  the  complex  amplitudes  Cp^(k)  satisfy  the  transversalify  and  reality  conditions 
which  are,  respectively, 

k-C?>(k)=0,  C<*>(k)  =  cf >(-k)  ■ 

The  representation  (12.1)  and  (12.2)  holds  for  any  moment  of  time.  Hence,  we  can  set  t  =  0 
to  generate  suitable  initial  conditions  for  an  electromagnetic  pulse  propagating  in  empty 
space  by  choosing  specific  functions  Co'^^(r). 

Consider  a  few  examples.  Let  there  be  translational  invariance  along  the  y  axis.  In  this 
case  the  fields  depend  orily  on  x  and  z,  i.e.,  r  =  {x,0,z).  Accordingly,  the  wave  vector  has 
the  form  k  =  {kx,  0,  kx)  and  dk  =  dkxdk^.  Let  62  =  (0, 1, 0)  be  the  unit  vector  along  the  y 
axis.  Introduce 

C±(r)  =  §62  Ao  gTiko-  =  62  ^^(r)  ,  (12.3) 

where  Aq  and  k  are  real  constants,  and  ko  is  a  fixed  wave  vector.  Note  that  the  field 
is  automatically  transversal.  The  corresponding  fields  determine  suitable  initial  conditions 
to  generate  a  pulse  propagating  in  the  direction  of  ko  whose  frequency  band  is  centered  at 
(Jq  =  cko  and  its  width  is  proportional  to  ck.  The  pulse  is  linearly  polarized  along  the  y  axis: 

Eo  =  ezFo,  Eo  =  C'^^-hC'(-^  (12.4) 

Bo  =  ^  V  X  62  ((7$+)  -  .  (12.5) 

The  action  of  the  diflferential  operator  is  defined  via  the  fast  Fourier  transform  (see  Section 
3)  on  the  grid  fine  enough  to  support  the  bandwidth  limited  function  (12.3).  In  the  Fourier 
basis,  iV/V— A  — >  —k/k. 
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To  obtain  suitable  initial  conditions  for  a  pulse  propagating  in  the  direction  ko  and  whose 
polarization  li^  in  the  xz-plane,  we  make  use  of  the  electromagnetic  duality  of  Maxwell  the¬ 
ory  which  states  that  the  dynamics  remains  unchanged  when  electric  and  magnetic  charges 
switch  places  and  simultaneously  Et  — >  — Bt  and  B*  —>■  Et.  According  to  the  duality  theo¬ 
rem,  we  can  take 

Eo  =  ^  Vxd2(C’{+)-C7S-))  (12.6) 

Bo  =  —©2  (^0^^  +  ^0  •  (12.7) 

Finally,  suitable  initial  conditions  for  a  pulse  propagating  in  the  direction  ko  with  a 
generic  elliptic  polarization  are  obtained  by  taking  a  linear  combination  of  the  above  two 
initial  conditions  for  two  independent  linear  polarizations  of  the  pulse.  The  amplitudes 
(r)  can  also  be  set  munericaUy  from  actual  measurements  of  a  particular  pulse  of  interest. 

II  Conductivity  and  absorbing  boundary  conditions 

In  numerical  simulations,  the  grid  in  coordinate  space  is  of  necessity  finite.  In  scattering 
problems  we  are  interested  in  the  pulse  shape  and  polarization  which  we  wish  to  compute 
in  the  asymptotically  large  coordinate  region.  This  requires  that  not  only  the  leading  edge 
of  the  reflected  pulse  should  have  reached  the  asymptotic  region,  but  also  the  trailing  edge 
should  have  done  so  as  well.  This  is  essential  if  the  reflected  pulse  propagates  in  a  highly 
dispersive  medium,  or  the  target  has  a  complex  shape,  or  both.  A  complication  arises  from 
the  very  nature  of  the  fast  Fourier  transform  method.  The  method  is  designed  to  describe 
periodic  functions  and,  consequently,  if  the  pulse  has  a  finite  amplitude  at  the  edge  of  the 
grid,  this  finite  value  would  appear  back  at  the  other  edge,  with  totally  disastrous  results 
for  the  computation.  In  quantum  computational  physics  this  problem  is  often  solved  by 
using  an  optical  potential  that  absorbs  the  signal  as  it  reaches  the  grid  boundary.  A  similar 
method  can  be  developed  for  our  treatment  of  Maxwell’s  theory.  Before  we  do  so  let  us 
point  out  that  an  absorbing  boundary  condition  is  not  the  only  way  to  solve  the  problem. 
For  instance,  in  the  case  of  a  complex  target,  an  aneillary  grid  may  be  defined  in  one  of  the 
coordinates  which  extends  to  large  distances.  The  pulse  may  be  transferred  in  a  gradual 
manner  firom  the  small  grid  (near  the  target)  to  this  larger  grid  to  prevent  the  pulse  from 
ever  reaching  the  edge  of  the  small  grid.  This  technique  can  also  be  applied  to  generate  a 
pulse  by  an  antenna  of  a  complex  construction.  The  dynamics  of  the  portion  of  the  pulse  on 
the  larger  grid  may  be  treated  analytically  (if  dispersion  properties  of  the  medimn  are  not 
too  complex). 

In  quantum  mechanics  absorbing  boxmdaxy  conditions  are  made  by  adding  an  imaginary 
potential  to  the  Hamiltonian  with  support  near  the  grid  edges.  In  the  Maxwell  theory,  the 
same  can  be  achieved  by  adding  conductivity  which  gradually  increases  as  the  grid  edges  are 


40 


approached.  An  interaction  of  conducting  media  with  electromagnetic  radiation  is  described 
by  Ohm’s  law, 

Jf  =  ,  (12.8) 

combined  with  Maxwell’s  equation  (2.1),  where  the  displacement  current  is  amended  as 
dtDt  +  (47rcr/c)Et  with  a  =  (r(r)  being  the  conductivity  of  the  medium.  Consider  a 
linearly  polarized  plane  wave  moving  along  the  z  axis.  Let  be  the  Fourier  transform  of  the 
only  component  of  the  electric  field  Et.  Disregarding  for  a  moment  any  possible  anomalous 
dispersion  of  the  medium,  we  find  that  satisfies  the  equation 


E^{z)  =  0  . 


(12.9) 


Equation  (12.9)  is  identical  to  the  stationary  Schroedinger  equation  with  an  optical  (absorb¬ 
ing)  potential  being  proportional  to  0{z).  In  simulations  of  quantum  wave  packets  it  has 
been  found  that  one  of  the  optimal  potentials  has  the  form  [17] 

(r(z)  =  {n  +  1)-'  (T„  (z/L)^  ,  n  >  2  ,  (12.10) 

in  the  interval  z  €  [0,  L]  and  o{z)  —  0  otherwise.  So,  our  next  task  is  to  find  an  optimal 
constant  cr„  such  that  that  the  conducting  layer  would  not  reflect  or  transmit  electromagnetic 
energy  in  some  designated  frequency  band.  Maxwell’s  equations  in  a  conducting  medium 
are  form-invariant  under  the  scaling  transformations 

a;  -*■  ,  <T(r)  ^u(/3r)  ,  ^J(r)  -»•  ,  (12.11) 


where  a  and  yS  are  positive  constants.  If  the  conductivity  <T(r)  was  found  optimal  for  a 
firequency  u  and  over  a  length  L,  then  the  optimal  conductivity  for  a  firequency  would 
be  0(r{0T)  and  the  new  length  over  which  it  is  taken  to  act  would  be  I-//3. 

Let  us  first  study  the  reflectivity  of  the  absorbing  layer.  Suppose,  (t{z)  =  ffod{z)  where 
6(z)  is  the  Heaviside  function.  If  a  monochromatic  linearly  polarized  wave  coming  from  the 
negative  z  region  has  an  amplitude  one,  then  the  reflected  wave  has  the  amplitude  [19] 


1 

l  +  i'u,  ’ 


The  energy  of  the  reflected  wave  is 


1/^  =  1  + 


4m<7Q 

UJ 


=  1  +  • 


i?a>  =  l?7a»P  =  gS/4  +  0(a;^)  . 


(12.12) 


(12.13) 


Ru  increases  as  the  ratio  gets  higher  and  is  small  if  9^/4  1.  Consider  now  (t{z)  which 

monotonically  increases  from  ^  =  0  in  the  positive  z  direction.  Let  ku  —  ku{z)  be  a  local 
wave  vector  of  the  wave  in  the  conducting  medium, 

k^{z)  =  c"^a?>/l  +  iqo,{z)  .  (12.14) 
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We  shall  argue  that  the  reflection  is  negligibly  small  if  the  local  wave  vector  does  not  signif¬ 
icantly  changes  over  a  distance  of  order  that  is, 


k^j{z  -h  Sz)  —  kujz) 

kM) 


«1 , 


Sz  ~  - 


(12.15) 


Making  a  linear  approximation  in  (12.15),  we  infer  that 

|5i9w(^')|  (2c)“^a>  (l  +  g^W)  ^  ’  (12.16) 


which  must  be  valid  for  all  values  of  qu,(z)  including  small  ones  when  the  reflection  is 
small.  Inequality  (12.16)  allows  us  to  reverse  the  argument,  that  is,  the  reflection  is  small 
if  \dzq^{z)\  «i  (2c)"^w.  Let  Oi  be  an  average  conductivity  over  a  layer  of  width  L, 
cri  =  L~^  dz(T(z).  In  particular,  for  (12.10),  -  (r„.  For  a  monotonically  increas¬ 

ing  function,  the  derivative  can  be  approximated  as  |529w(2;)|  ^  o'i/L.  This  leads  to  a 
necessary  condition  on  conductivity  to  suppress  the  reflection,  namely. 


o-L  < 


uPL 

8tc 


(12.17) 


Our  analysis  is  valid  if  the  higher  derivatives  of  {t(z)  are  not  large.  This  condition  requires 
that  the  exponent  n  in  (12.10)  should  not  be  less  than  two  to  insure  a  smooth  behavior  at 
2  =  0. 


The  transmission  can  be  estimated  as  follows.  Suppose  the  pulse  occupies  a  compact 
region  Ci.  Let  £p  be  the  pulse  energy.  The  pulse  looses  its  energy  as  it  propagates  through 
a  conducting  medimn  according  to  Ohm’s  law,  so  that 

c-%e^  =  -(2c)-'  f  dr  (tE^<  -87r(Tn5f  ,  (12.18) 

Ju 

where  an  =  niaxn  a.  Therefore  the  pulse  energy  decay  can  be  bounded  from  above  by 


<  e-8^n  ^ 


(12.19) 


In  the  one  dimensional  case  (12.10),  an  =  ai/{n  +  1).  For  the  time  t  —  Lfc  needed 
for  a  pulse  to  get  through  the  layer  of  width  L,  the  attenuation  should  be  large,  that  is, 
8T^Lahfc{n  -f  1)  ?$>  1.  Thus,  the  necessary  conditions  to  suppress  both  transmission  and 
reflection  (that  is,  to  ensure  an  almost  total  absorption)  of  the  pulse  arc 


{n  +  l)c 

SirL  ^  Sttc 


(12.20) 


By  changing  the  Hamiltonian  yP,  the  conducting  layer  can  be  included  into  the  split 
or  leapfrog  algorithm.  Since  the  conducting  layer  produces  attenuation,  the  conductivity  a 
must  be  included  into  the  operator  C^t  hi  the  modified  leapfrog  scheme.  It  is  also  possible 
to  create  an  absorbing  and  non-reflecting  layer  by  using  a  passive  medium  (e.g.  a  Lorentz 
model).  The  analysis  of  the  medium  properties  would  be  similar  to  that  for  a  conducting 
layer.  In  fact,  using  a  layer  of  a  passive  medium  would  offer  more  flexibility  in  solving  the 
grid  boundary  problem. 
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l lb.  Enter  the  indirect  cost  rate  in  effect  during  the  reporting 
period. 

llc.  Enter  the  amount  of  the  base  against  which  the  rate 
was  applied. 

l ld.  Enter  the  total  amount  of  indirect  costs  (^rged  during 
the  report  period. 

l le.  Enter  the  Federal  share  of  the  amount  in  1 1d. 

Note:  If  more  than  one  rate  was  in  effect  during  the  period 
shown  in  item  8,  attach  a  schedule  showing  the  bases 
against  which  the  different  rates  were  applied,  the 
respective  rates,  the  calendar  periods  teey  were  in 
effect,  amounts  of  indirect  expense  charged  to  the 
project,  and  the  Federal  share  of  indirect  expense 
charged  to  the  project  to  dale. 
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Abstract 

Under  rather  general  assumptions,  and  in  a  relatively  simple  and 
streughtforward  manner,  it  is  shown  that  the  characteristics  of  signals 
which  travel  through  homogeneous,  as  well  as  inhomogenous,  passive 
media  have  the  principal  features  usually  associated  with  the  phenom¬ 
ena  of  precursors,  as  generally  follows  from  more  elaborate  studies. 
The  simplicity  of  the  present  arguments  permit  analytic  studies  to  be 
made  for  a  greater  variety  of  media  than  is  normally  the  case. 


Introduction 

The  propagation  of  electromagnetic  signals  through  passive  media,  i.e.,  media 
exhibiting  both  dispersive  and  absorptive  properties,  of  many  different  kinds 
is  a  problem  of  interest  in  many  applications.  For  simplicity  of  analysis,  we 
model  such  problems  in  relatively  economical  but  fairly  general  terms.  Let 
fo{t)  denote  the  temporal  dependence  of  the  initial  signal  at  some  initial 
depth  level  z  0  inside,  or  at  the  surface,  of  some  medium.  We  presume  the 
medium  to  be  homogeneous  at  this  stage  and  that  its  response  characteristics 
are  stationary.  We  denote  by  fz{t)  the  temporal  dependence  of  the  signal 
after  traveling  a  distance  z  >  0  through  the  passive  medium.  The  relation 
between  the  input  and  output  signals  is  that  of  a  standard  filtration  which 
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may  be  presented  in  several  diflPerent  forms.  One  such  form  is  given  by 
convolution  with  the  medium  impulse  response  function  mz{t),  and  takes 
the  form 

fz{t)  =  j  s)fo{s)  ds  .  (1) 

Note  that  all  integrals  with  unspecified  limits  are  assumed  to  extend  over 
the  whole  real  line.  As  a  convolution,  these  functions  have  a  multiplicative 
relation  in  frequency  space  given  by 

F,{u;)  =  M,{u;)Fo{uj),  (2) 

where 

Fz{uj)  =  Je^Mt)dt,  (3) 

Mz{uS)  =  J  e^mz(t)dt.  (4) 

As  a  causal  phenomena  (no  signal  out  before  a  signal  in),  it  follows  that 
Tnz{t)  =  0,  <  <  0,  a  property  that  means  that  Mz{(jS)  has  no  poles  in  the 
upper  half  complex  uj  plane,  i.e.,  no  poles  for  a;  =  a;r  +  i^Oi  with  Ui  >  0. 
These  few  properties  are  standard  for  causal  signal  transmission  [1]. 

Let  us  now  take  up  an  elementary  property  of  Mz{oj)  regarding  its  de¬ 
pendence  on  the  depth  parameter  z.  Clearly, 

{^) 

=  MzA^)MzMFo{oj) 

—  Mzi+z2i^)Fo{u>) .  (5) 

Since  Fo{uj)  is  arbitrary  in  this  expression,  it  follows  that 

Mzi+z2{'^)  =  Mzi{u})Mz2{i^)  (6) 

holds  for  all  nonnegative  Zi  and  Z2.  In  particular,  it  follows  that 

Mo(aj)  =  Mo{ui)^  (7) 

implying  that  either  Mo(u})  =  1  or  Mo{oj)  =  0.  If  there  are  frequency  bands 
that  Fo(a;)  will  never  jKSsess,  then  it  is  possible  to  choose  Mo{uj)  =  0  in  that 
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band,  but  for  simplicity  we  may  as  well  assume  that  Mo{uj)  =  1  for  all  u 
since  the  final  answer  will  be  imchanged.  Moreover,  on  physical  grounds  it 
is  natural  to  assume,  for  each  a;,  that  iV4(ct?)  is  a  continuous  function  of  2:. 
In  that  case,  a  standard  argument  asserts  that  the  only  function  that  fulfills 
(6)  necessarily  has  the  functional  form 

MM  = 

Let  us  recall  that  argument. 

To  establish  the  stated  form  for  we  first  observe  that 

MM  =  [M.iM]"  (9) 

for  an  arbitrary  positive  integer  n.  Since  Mz{ui)  is  assumed  to  be  continuous 
and  Mo(a7)  =  1,  it  follows  that  MzM)  ^  0  for  large  enough  n,  and  so 
Mz{oS)  ^  0  for  any  2.  Next,  let  Rz{(jJ)  =  ln(M^(a;)),  which  is  well  defined 
for  all  z  and  a;;  furthermore,  (6)  leads  to  =  Rzii^)  +  Rz2i^)-  This 

relation  implies  that  Rnziu)  =  nRz{oj)  as  well  as  Rz/m{^)  —  {\/f^Rz{^), 
for  arbitrary  positive  integers  n  and  m.  Consequently, 

Rnz/mioj)  =  {nJm)Rz{oS) .  (10) 

Now  let  =  1  so  that  Rn/M)  =  {n/m)Ri{u))\  then  take  a  limit  such  that 
(n/m)  z,  an  arbitrary  positive  number.  By  the  assumption  of  continuity, 
this  limit  converges  and  we  learn  that 

Rz{u))  =  zRi{uj)  =  -zA(u!) .  (11) 

This  concludes  the  proof  of  (8). 

As  the  transfer  function  of  a  passive  medium  it  is  necessary  that  |  Mz  {ui)  |  < 
1  in  which  case  we  also  find  that  ReA(a;)  >  0.  To  simphfy  notation,  we  let 

A{oj)  =  Biu;)  +  iC{u),  (12) 

where  B{oj)  =  ReA(a;)  and  C{(jS)  =  ImA(c<;).  Thus  we  require  that 

B{u;)  >  0  .  (13) 

By  causality,  B{(jj)  and  C{u)  are  related  by  Hilbert  transforms,  while  C{uj)  is 
undetermined  by  this  procedure  up  to  the  contribution  of  a  so-called  Blaschke 
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factor.  We  also  note  that  for  real  functions  /^(t)  for  all  2:  and  any  initial  signal 
it  necessarily  follows  that  the  impulse  response  function  Tnz{t)  is  a  real 
fimction.  Consequently,  Mz{-u))  =  which  implies  that 

=  (14) 

C(-uj)  =  -C{uj)  .  (15) 

Since  B{u})  >  0  it  follows  that  the  function  B{ui)  has  special  properties  at 
and  near  any  frequency,  say  u/,  for  which  5  (a/)  ==  0.  Assuming  on  physical 
grormds  that  B(uj)  is  continuous  in  u),  and  even  has  a  few  derivatives  [e.g., 
B{(J)  G  C%  then  the  behavior  of  B  in  the  vicinity  of  u/  is  given  by 

B(u;)  =  ^a-^(u-u/y,  (16) 

where  a  >  0  has  the  imits  of  acceleration  (e.g.,  m/s^). 

Let  us  focus  on  the  behavior  of  A{uj)  in  the  vicinity  of  =  0.  In  that 
case,  and  to  sufficient  accuracy,  we  have 

B{u))  =  L~^  +  ,  (17) 

C{uj)  = —v~^u; ,  (18) 

where  L,  v,  and  a  denote  constants  with  the  dimensions  of  len^h,  velocity, 
and  acceleration,  respectively.  The  behavior  at  w  =  0  forces  >  0.  More¬ 
over  when  S(0)  =  L~^  =  0,  it  follows  that  a~^  >  0;  in  that  case,  we  shall 
adopt  the  generic  form  and  assume  that  a  >  0. 

It  is  natural  to  ask  the  question:  Under  what  circumstances  does  it  follow 
that  B{0)  =  0?  A  common  model  of  a  passive  medium  is  based  on  the 
Lorentz  model  [2].  The  physical  principle  behind  this  model  is  the  response 
of  a  damped  harmonic  oscillator  to  a  driving  force.  The  physical  model  in 
this  case  is  given  by 

mq{t)  +  bq{t)  +  kq(t)  =  h{t)  ,  (19) 

where  m  >  0,  6  >  0,  and  fc  >  0  denote  the  inertial,  damping,  and  restora¬ 

tive  parameters,  q{t)  denotes  the  oscillator  amplitude,  and  h(t)  denotes 
the  applied  force.  Assuming  that  any  transients  have  decayed  away,  let 
us  focus  on  the  steady  state  solution.  If  h{t)  =  U  cos{u)t),  then  q{t)  = 
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V  cos(a;i)  +  W  siii(a;<),  where 


(fc  —  map)  U 

{k  —  maPY  +  ' 

bivU 

(fc  —  mupy  +  ipaS)^ 


(20) 

(21) 


If  the  driving  frequency  u)  =  0  for  this  model  it  follows  that  V  —  U/k  and 
W  =  0.  Since  this  result  is  entirely  independent  of  the  dissipation  coefficient 
6,  it  holds  even  when  6  =  0,  namely,  it  holds  in  the  case  there  is  no  loss 
whatsoever.  Consequently,  for  such  a  model  it  follows  that  at  a;  =  0  there 
can  be  no  loss  and  so  B{0)  =  0.  Of  course,  such  an  argument  does  not  apply 
to  a  metal,  but  it  does  hold  for  a  broad  spectrum  of  materials. 

Although  we  have  argued  that  B{0)  =  0  for  a  specialized  model,  it  is 
quite  plausible  that  such  a  result  holds  for  a  wide  class  of  models  as  well. 
One  such  more  general  model  would  be  of  the  form 


+  6(9, (()  +  S/.1  9j(*)  =  fc(i) ,  l<i<J,  (22) 

which  corresponds  to  a  situation  where  h{t)  acts  as  the  driving  force  on  a 
collection  of  coupled  damped  harmonic  oscillators.  Again  we  may  consider 
the  steady  state  response  when  h(t)  =  U  cos(a;t).  In  particular,  if  o;  =  0  the 
solution  is  totally  independent  of  any  dissipation  parameters.  Consequently, 
it  follows  that  B(0)  =  0  in  this  case  as  well.  It  is  not  difficult  to  imagine  other, 
more  complicated  situations  for  which  the  same  result  holds.  Therefore,  on 
the  basis  of  a  wide  class  of  physically  plausible  models,  we  are  encouraged 
to  consider  the  consequences  of  the  assumption  that  B{aS)  =  ^a~^uP,  o  >  0, 
when  oj  is  small. 

We  now  examine  the  signal  behavior  as  a  function  of  z  that  follows  from 
the  foregoing  discussion.  In  particular,  we  first  focus  on  the  impulse  response 
function 


^ziv-za  ^w^l2-zA(w)  ^ 


(23) 


where  A{uj)  =  B{aS)  +  iC{aj),  and  B{aS)  =  0(aJ^)  while  C{a})  =  0{uP).  To 
simplify  matters,  let  us  further  assume  that  B{aj)  is  boimded  away  from  0  for 
nonzero  values  of  cj.  In  that  case,  for  large  z,  we  observe  that  the  integrand 
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tends  to  be  strongly  suppressed  for  u;  away  from  zero.  To  sufficient  accuracy 
it  follows  that  we  may  set 


Consequently,  for  large  z  values  we  learn  that 

fzit)  =  J  mz{t-  s)fo{s)  ds 

I  a  f  (  a{t-  s-  v~'^zy 


There  are  two  characteristic  features  of  this  resultant  behavior  that  are 
worth  emphasizing: 

1.  According  to  the  square  root  prefactor,  it  follows  that  the  amplitude 
decays  proportional  to 

2.  According  to  the  integral  factor,  the  result  has  a  temporal  depen¬ 
dence  determined  by  the  form  of  the  original  signal  centered  at  a  temporal 
value  oft-v~^z,  and  then  integrated  over  that  region  weighted  by  a  Gaussian 
of  increasing  width. 


Approximate  causality 

Before  passing  to  several  examples  of  (25)  it  is  important  to  observe  that  the 
medium  impulse  response  function  (24)  is  not  causal  since  it  does  not  vanish 
for  t  <  0.  In  point  of  fact,  the  expressions  (17)  and  (18)  do  not  make  a 
strictly  causal  impulse  function.  However,  they  make  an  approximate  causal 
function  provided  that  az  »  since  in  that  case  the  impulse  function  is 
extremely  small  for  t<0.  Let  us  examine  this  situation  more  closely. 

In  order  for  the  medium  function  A{cj)  to  correspond  to  a  causal  function 
it  shoiild  be  represented  in  the  form 

fOO 

A{uj)=  /  [l-e“^^]fc(s)ds,  (26) 

Jo 
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where  >1(0)  =  0,  and  in  order  that  Re^(w)  =  B{uj)  >  0  it  is  sufficient 
that  the  weight  function  k{s)  >  0.  Observe,  by  this  representation,  that 
A{u})  is  analytic  in  the  upper  half  plane,  which  is  a  necessary  and  sufficient 
condition  for  the  impulse  function  to  be  causal.  Now,  as  an  example,  consider 
k{s)  =  where  K'  and  K  are  positive  parameters  to  be  determined. 

In  that  case 


poo 

j5(a>)  =  K'  [1  -  cos(a>s)]  ds  , 

Jo 

(27) 

poo 

=  1  sm{ujs)e~^*  ds  . 

Jo 

(28) 

We  are  intested  in  a  large  value  for  both  K  and  K'  such  that,  to  some 
approximation,  k(v)  acts  rather  like  a  delta  function  in  the  integrand  focusing 
on  values  of  v  near  zero.  Thus,  to  a  certain  approximation. 

B{uj)  =  s'^e  - , 

(29) 

i*CO 

c{uj)  -  -K'u}  /  ds  +  •  •  •  . 

(30) 

On  comparison  with  (17)  and  (18)  we  see  that 

roo 

=  5^6”^"  ds  -  2K'/K^  , 

Jo 

(31) 

Ts  ds  -  K'/K^  . 

Jo 

(32) 

Thus  a  =  \Kv.  Provided  >  1,  this  approximation  is  satisfactory  from 
a  practical  viewpoint;  otherwise  some  of  the  additional  terms  in  (29)  and 
(30)  need  to  be  taken  into  account.  Unfortunately,  introducing  additional 
terms  makes  the  medium  impulse  response  function  no  longer  anal3dically 
tractable,  generally  speaking. 

Specific  examples 

Let  us  illustrate  the  resultant  signal  (25)  in  two  important  examples.  First  we 
consider  an  analytically  simple  case,  namely  a  Gaussiarifinitikl  pulse.  Thus, 
let 

fo{t)  =  cos{ujot) ,  (33) 
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where  we  refer  to  T  as  the  pulse  width.  In  this  case,  (25)  reads 


g  a(s  t+v  ^z)'^/2z^  a^l2T'^  COs{u}oS)  ds  . 


An  elementary  integration  leads  to 


2(1  + af^ - V(l  +  ./aTJ-<^> 

If  we  assume  that  2  ^  aT^,  then  in  effect  (35)  reduces  to 


a(t  —  v  ^zY  {i^pTy 


The  peak  amplitude  of  the  resultant  pidse  /^(t)  occurs  when  z  =  vt,  or  stated 
otherwise  when  t  —  v~^z.  At  that  point  in  z  (or  at  that  moment  in  time 
t),  the  amplitude  has  been  reduced  from  its  initial  maximum  of  unity  by  a 
factor  proportional  to  z~^l'^.  The  resultant  pulse  shape  in  the  present  case 
is  still  Gaussian,  but  it  is  substantially  broadened  since  the  pulse  width  T 
of  the  input  Gaussian  has  been  replaced  by  ^Jzja  T.  If  we  assume  the 
energy  (or  power)  in  a  given  pulse  fz{t)  is  proportional  to  J\fz{t)\^  dt,  then 
it  follows,  for  z  »•  aT^,  that  the  normalized  signal  energy  is  given  by 

(37) 

/i/o(t)Pdt  y  z 

As  anticipated,  this  ratio  decreases  as  z  increases,  but  thanks  to  the  pulse 
broadening,  it  decreases  only  as  z~^^^  rather  than  z~^,  as  may  have  been 
expected. 

As  a  second  example,  we  suppose  that 

/o(f)  =  rect(t/T)  coe{uJot)  ,  (38) 

which  leads  to  the  output  signal 

Mt)  =  ^-a{s-t+v-^zfl2z  cos(a;„s)  ds  .  (39) 

y  I-Kz  y_T/2 
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This  integral  depends  on  z  in  a  moderately  involved  way,  and  for  convenience 
we  only  evaluate  (39)  in  the  case  that  z  »  aT^.  Then  we  find 


=/ 


~a{t—v  '^zYl2z 


cos,{uJos)  ds 


aT^  rsin(a?or/2)l 
27rz  L  i^oT/'^  J 


(40) 


Again,  the  result  is  that  of  a  traveling  pulse  moving  essentially  at  the  speed 
V,  the  amplitude  of  which  falls  oflF  with  distance  like  z  provided  of  course 
that  sin(a;oT/2)  ^  0.  For  the  case  that  sin(a;oT/2)  =  0,  iVoT  ^  0,  see  below. 


More  general  inputs 

Apart  from  an  unimportant  factor,  observe  that  the  amplitude  factor  from 
/o(t)  in  (36),  i.e.,  exp[-a;^TV2],  is  basically  Fo(0)  -  / /o(s)  ofe, .  and  the 
amplitude  factor  in  (40),  i.e.,  sin(a;oT/2)/(a;oT'/2)  is  essentially  just  Fo(0)  = 
J  fo{s)ds.  We  can  understand  and  generalize  this  result  as  follows.  Since 

fz{t)  =  Jm,it-s)Ms)ds,  (41) 

it  follows,  for  an  initial  signal  which  is  much  narrower  in  the  temporal  do¬ 
main  than  the  medium  response  function,  that  we  can  make  a  Taylor  series 
expansion  of  the  medium  response  function.  In  particular,  we  have 

=  J['>T^z{t)fois)  -  m'^{t)sfo{s)  +  |m"(t)s2/o(s)  -k  •  •  •  ]  ds 
=  m^(<) / fo{s)ds  -  m!^{t)Jsfo{s)ds  +  |m"(t)/sVo(s)ds  +  •  •  •  . 

(42) 

When  m^{t)  is  sufficiently  broad  in  time,  then  its  time  derivatives  are  typ¬ 
ically  much  smaller,  and  the  indicated  sum  involves  terms  that  are  rapidly 
diminishing  relative  to  the  former  terms. 

As  an  example  of  the  application  of  (42),  suppose  that  tVoT  =  2n7r,  n  € 
{1, 2, 3, . . then  (40)  vanishes  because  ffo(s)ds  =  0.  Moreover,  since  fo(t) 
in  (38)  is  an  even  function,  it  also  follows  that  f  s  fo{s)  ds  =  0.  Thus  according 
to  (42),  we  find  to  leading  order  that 

fz(t)  =  im':{t)Js^fo{s)ds,  (43) 
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namely, 


where  the  amplitude  has  been  evaluated  when  lJoT  =  2n7r. 


(44) 


Commentary 

Observe,  that  although  the  original  Gaussian  and  the  rectangular  envelopes 
are  fundamentally  different,  the  pulses  that  result  after  a  long  propagation 
distance  aT^  are  essentially  the  same,  differing  only  in  their  respective 
amplitudes.  Those  amplitudes  effectively  measure  the  spectral  content  of 
the  original  pulse  at  zero  frequency  where  the  narrow  window  of  transmis¬ 
sion  exists,  and  for  the  Gaussian  envelope  that  content  is  much  smaller  than 
is  available  with  the  rectangular  envelope.  This  result  is  entirely  consistent 
with  the  concept  that  the  output  selects  its  preferred  signal  shape  consistent 
with  the  frequency  window  allowed  by  the  medium.  Ultimately,  this  means 
that  the  transmitted  signal  is  Gaussian  in  character  with  an  ever  increasing 
pulse  width  in  accord  with  the  ever  narrowing  Gaussian  pass  band  as  z  in¬ 
creases.  This  kind  of  temporal  behavior  is  essentially  the  kind  of  behavior 
typically  ascribed  to  the  so-called  Brillouin  precursor  [2].  It  is  paradoxical, 
but  true,  that  Brillouin  precursors  require  pulse  envelopes  with  high  fre¬ 
quency  content.  Sharp  edges  in  pulse  envelopes  significantly  enhance  the 
high  frequency  content  as  compared  to  smooth  envelopes,  and  this  situation 
is  clearly  illustrated  in  Eqs.  (36)  and  (40). 

However,  sharp  edges  in  pulse  envelopes  are  not  the  only  way  to  raise  the 
high  frequency  content. 


Chirp  Signals 

As  an  initial  signal  consider  the  linearly  chirped  signal  [3] 

fo{t)  =  e{t)  coe{uJot  +  ,  (45) 

where  the  signal  envelope  e{t)  has  a  pulse  width  of  T  and  we  assume  that 
aT^  ^  1.  The  Fourier  transform  of  (45)  may  be  approximately  given  by 


(46) 


V  ot 


The  frequency  content  at  a;  =  0  is  then  given  by 


As  an  illustration,  consider  e(t)  =  exp(— <^/2r^).  In  this  case, 


Fo(0)  = 


•ia^T* 


l["-(S=) 


+  sin 


( 


aT^  )i  ' 


(47) 


(48) 


As  noted  above  this  factor  essentially  dictates  the  amplitude  for  the  output 
Gaussian  in  (24). 

The  main  point  to  observe  in  this  calculation  is  the  significant  enhance¬ 
ment  of  the  overall  amplitude  of  the  signal  that  will  emerge.  In  particular, 
the  primary  change  has  been  that  the  factor  exp[— (a;oT)^/2]  in  (36)  has  been 
replaced  by  exp[-{oJoTy where  aT^  >  1.  This  model  calculation 
strongly  suggests  that  combining  chirp  signals  with  envelopes  having  sharp  or 
nearly  sharp  edges  should  lead  to  signals  with  relatively  significant  precxusor 
content. 


Composite  Media 

Up  to  this  point  we  have  assumed  that  the  signal  propagates  through  a 
single,  homogenous  medium,  the  properties  of  which  are  characterized  by  the 
complex  medium  function  A{u}).  Suppose  instead  we  deal  with  a  composite 
media  made  up  from  J  layers  each  of  thickness  l\  =  z\,  l2  =  Z2  —  zi,...,lj  — 
zj  —  zj-i.  To  each  layer  we  attribute  a  medium  function  Aj{oj),  1  <j  <  J- 
If  we  assume  that  Zj  <  z  <  Zj+i,  for  example,  then 

Mz(u)  =  exp[-(2  -  Zj)Aj+i{uj)  -  hAkiuj)] .  (49) 

A  more  convenient  expression  for  such  media  is  given  by 

Mz{u})  =  exp[—JlAu{u!)du]  ,  (50) 

which  applies  for  piecewise  constant  functions  A„(ct;).  It  is  plausible  that 
such  an  expression  applies  as  well  for  slowly  varying  media  functions  A„(a;). 
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As  was  previoTJsly  the  case,  it  is  iiseful  to  separate  A„(a;)  into  real  and 
imaginary  parts,  i.e., 

Mz{ui)  —  exp[— du  —  iflCu{(*^)  du] .  (51) 

Under  suitable  conditions  we  may  set 

Bu(uj)  =  ,  Cu{u^)  =  -v-'^{u)uj  ,  (52) 

in  which  case  we  are  led  to 

Mz{u))  =  exp[— |/oa“^(«)a>^ du  +  i J^v~^{u)u;du]  .  (53) 

Equations  (51)  and  (53)  apply  to  various  situations. 

As  a  simple  example  of  how  these  equations  may  be  applied  let  us  con¬ 
sider  a  signal  that  penetrates  a  passive  medium  with  a  finite  thickness  of  £ 
and  otherwise  propagates  in  free  space.  Therefore,  in  this  case  >  0 

for  0  <  w  <  ^,  while  for  u  >  £,  Bu{oj)  =  0.  This  example  pertains  to  a 
radar  signal  that  penetrates  a  canopy  of  foliage  or  some  other  absorbing  and 
dispersive  medium.  Even  a  comparatively  thin  passive  medium  may  be  well 
approximated  by  (53)  when  A„(a;)  is  sufficiently  large  for  a;  0.  Under  the 
stated  conditions  it  follows,  when  z>  £,  that 

Mj(u;)  =  exp{— ^^/2a  +  i[{z  —  £)ul c  +  £u)/v])  ,  (54) 

where  c  is  the  speed  of  light  in  vacuum,  and 

d=£ [/^^(tt)”^ d'u\~^  , 

v  =  e[^[v{u)-Uu]-K  (55) 

Since  v(u)  <  c,  it  follows  that  v  <  c.  For  simplicity  in  our  example,  let  us 
assume  that  v  c.  Finally,  when  £/a  is  sufficiently  large,  we  are  led  to 

fz{t)  =  V^£  JMs)  ds  ,  (56) 

defined  for  z  >  £.  If  Fo(0)  =  0,  then  appeal  can  be  made  to  (42). 

This  equation  should  be  compared  with  Elq.  (36)  that  corresponds  to 
propagation  within  a  single,  homogeneous  passive  medium.  The  temporal 
signal  in  both  cases  is  Gaussian  and  apart  from  the  overall  amplitude,  the 
shape  is  principally  determined  by  the  passive  medium.  This  result  holds 
whenever  the  medium  acts  as  a  low  pass  filter  for  an  initial  signal  that, 
comparatively  speaking,  has  a  broad  spectral  content  including  a;  =  0. 
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Stochastic  Media 


Under  certain  cimcumstances,  it  is  appropriate  to  regard  the  passive  medium 
as  a  stochastic  variable  and  to  define  the  output  signal  as  the  ensemble  av¬ 
erage  of  a  stochastic  output.  A  description  of  this  kind  may  arise  if  the 
output  signal  is  the  superposition  of  multiple  signals  from  a  rapidly  chang¬ 
ing  collection  of  media  profiles.  More  specifically,  such  a  situation  may  arise 
in  the  case  of  a  coherent  superposition  of  multiple  radar  pulse  returns  that 
have  traversed  a  rapidly  fluctuating  medium.  While  details  of  the  fluctuating 
medium  may  be  too  diflBcult  to  describe,  certain  overall  features  —  stochas¬ 
tic  invariants  -  may  provide  gross  characteristics  with  which  selected  signal 
properties  may  be  analyzed. 

Let  us  focus  on  propagation  through  a  stochastic  but  otherwise  homoge¬ 
neous  medium.  In  mathematical  terms  this  means  that  Mz{uj)  and  thereby 
A{ijj)  are  regarded  as  suitable  stochastic  variables.  The  observed  signal  f°{t) 
is  determined  by  the  expression 

where  {(•))  denotes  an  ensemble  average.  Let  us  illustrate  the  kind  of  be¬ 
havior  that  may  occur.  For  that  purpose  we  again  suppose  that 

A{u)  =  (58) 

and  we  assume  that  a  is  the  only  random  variable.  As  an  example  we  suppose 
that  the  variable  o,  which  can  only  be  positive,  is  distributed  according  to  a 
probability  density  given  by 

p(a-i)  =  ^  ,  0  <  0"^  <  oo  ,  (59) 

where  6  is  a  parameter  which  essentially  sets  the  scale.  In  particular,  the 
mean  value  of  a~^  is  given  by 

p{a~^) da~^  —  -  .  (60) 

As  a  consequence,  we  find  that 


it-z/v)  ^oi^) 

(1  -I-  2;a;2/6)'”+i 


duj  . 


(61) 
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Even  without  evaluating  the  Fotuier  transform  it  is  clear  that  our  stochas¬ 
tic  example  has  led  to  a  dramaticsdly  enhanced  pass  band  at  low  firequencies. 
This  result  is  due,  of  course,  to  the  nonvanishing  probability  that  arbitrarily 
small  values  of  oT^  are  present  in  the  ensemble,  despite  the  fact  that  the 
probability  that  ==  0  is  strictly  zero  for  any  m  >  1. 

As  outlined  in  (42),  the  leading  behavior  of  (61)  is  given  by 

f°{t)  =  m%t)  Jfo{s)  ds  ,  (62) 


where 


m 


du) . 


(1  +  zuj^/b)"^+^ 

This  expression  may  be  evaluated  in  the  following  form: 


(63) 


m: 


=  (v^ I* -  ,  (64) 


where  the  coefficients,  for  0  <  1  <  m  and  1  <  m,  are  determined  by  the 
recursion  relation 

CM  =  I 

=  (2m  -  1)!! 

=  (2m  -  1  -  0  -h  1  <  1  <  m  -  1  .  (65) 


Conclusions 

Traditional  radar  applications  involve  signal  generation,  transmission,  and 
reception,  with  little  or  no  spectral  modification  save  for  a  possible  spec¬ 
tral  shift  (or  dilation)  due  to  relative  motion.  By  contrast,  radar  signals, 
or  electromagnetic  waves  more  generally,  that  traverse  passive  media  with 
normal  absorptive  and  dispersive  properties,  typically  suffer  significant  spec¬ 
tral  distortion.  For  example,  over  a  large  spectral  domain,  an  exponential 
attenuation  effectively  removes  any  signal  energy.  An  exception  to  exp>o- 
nential  attenuation  occurs  at  zero  frequency  when  attenuation  is  frequently 
absent.  On  the  basis  of  generally  applicable  arguments,  and  imder  a  wide 
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set  of  circumstances,  it  has  been  clearly  demonstrated  in  this  paper  that 
the  amplitude  of  the  transmitted  signal  decays  as  the  inverse  square  root  of 
the  distance  traveled  rather  than  exponetially.  Such  behavior  is  exactly  that 
characteristic  of  so-called  Brillouin  percursors  [1,  2].  While  such  phenomena 
are  jfrequently  discussed  in  terms  of  more  advanced  mathematical  techniques, 
it  is  demonstrated  in  the  present  work  that  they  can  be  largely  imderstood 
in  terms  of  fairly  simple  pass  band  arguments.^ 

The  simplicity  of  the  present  discussion  has  enabled  us  to  determine  the 
form  of  agnals  that  have  traversed  inhomogeneous  or  stochastic  media  as 
well.  The  robustness  of  the  characteristic  precursor  decay  rate  has  also  been 
established  for  these  more  general  media. 
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Abstract 

Maxwell’s  theory  for  an  electromagnetic  pulse  propagating  in  passive  linear  media 
is  reformulated  in  the  Hamiltonian  formalism  so  that  the  time  evolution  is  governed  by 
the  first-order  equation,  d^jdt  =  Wi,  where  'His  a  linear  differential  operator  (Hamil¬ 
tonian)  acting  on  a  multi-dimensional  vector  W  built  of  the  electromagnetic  inductions 
and  auxiliary  matter  fields  describing  the  medium  response.  The  initial  value  problem 
is  then  solved  by  means  of  the  time  leapfrog  method  in  combination  with  the  Fourier 
pseudospectral  method,  leading  to  a  time  domain  algorithm  for  numerical  simulations 
of  wide-band  electromagnetic  wave  packet  propagation  and  scattering.  The  algorithm  is 
applied  to  scattering  of  a  femtosecond  laser  pulse  on  resonant  transmission  and  reflection 
gratings  made  of  dispersive  (Drude  metals)  and  dielectric  materials. 
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1  Lorentz  models  in  the  Hamiltonian  formalism 

Propagation  of  an  electromagnetic  wave  packet  in  passive  (dispersive  and  absorbing)  media  in 
the  absence  of  external  radiating  sources  is  described  by  the  following  set  of  equations  [1] 

dtDt  =  cVxHt,  9tBt  =  -cVxEt,  (1-1) 

Dt  =  Et  +  Pt ,  Bt  =  Ht  +  Mt ,  (1*2) 

V-Dt  =  V-Bt  =  0,  {1-3) 

where  E*  and  H*  are  electric  and  magnetic  fields,  respectively,  D*  and  are  the  corresponding 
inductions,  Pt  and  Mt  are  the  medium  polarization  and  magnetization  vectors,  c  is  the  speed 
of  light,  dt  stands  for  the  partial  derivative  with  respect  to  time  t,  boldface  letters  denotes 
three-vector  fields  in  ^  whose  spatial  argument,  denoted  below  by  r,  is  suppressed,  and  the 
time  dependence  is  indicated  by  a  subscript.  The  time  evolution  of  the  medium  response  is 

usually  described  in  the  framework  of  a  particular  microscopic  medium  model.  In  the  present 

study,  multiresonant  Lorentz  models  are  considered.  A  general  case  within  the  formalism 
being  developed  can  be  found  in  [2].  In  the  Lorentz  model,  the  medium  magnetization  is  zero, 
while  the  medium  polarization  is  described  by  a  set  of  decoupled  second-order  differential 
equations  [1] 

N  ■ 

a?P“  +  2jadtPt  +  =  iv%Et ,  Pt  =  X) 

0=1 

where  Ua  s-re  resonant  frequencies,  ja  the  damping  coefficients,  and  Wpo  are  the  so  called 
plasma  frequencies. 

In  our  approach  no  boundary  conditions  are  imposed  on  electromagnetic  fields  at  medium 
and/or  target  interfaces.  The  latter  are  modeled  by  spatially  dependent  couplings  of  media 
to  electromagnetic  fields,  ujpa  =  Wpo(r).  At  any  interface,  the  couplings  are  allowed  to  have 
discontinuities,  or,  from  the  physical  point  of  view,  they  remain  smooth  but  change  rapidly, 
Au/iVcjpl/Wp  >?>  1,  at  the  interface,  where  A,„  is  a  typical  wave  length  of  the  incoming  wave 
packet.  The  conventional  boundary  conditions  are  automatically  generated  by  the  dynamics 
[!]• 

The  initial  value  problem  is  solved  in  the  space  of  square  integrable  functions  for  every 
component  of  the  medium  polarization  and  electromagnetic  fields.  This  implies  that  the  energy 
of  the  propagating  wave  padcet  remains  finite  (in  contrast  to  the  scattering  matrix  approach 
based  on  plane  wave  solutions).  The  initial  conditions  imposed  are  Pt=o  =  ^tP*=o  =  0 
the  electromagnetic  fields  form  an  initial  wave  packet  built  of  radiation  electromagnetic  fields. 

The  best  way  to  solve  the  initial  value  problem  for  a  system  of  differential  equations  is  to 
find  its  fundamental  solution.  The  latter  problem  is  often  studied  by  means  of  the  Hamiltonian 
formalism  [3].  In  the  framework  of  the  Hamiltonian  formalism,  an  original  system  of  differential 
equations  is  transformed  to  an  equivalent  system  of  first-oder  (in  time)  differential  equations  by 
expanding  the  original  configuration  space,  that  is,  by  going  over  to  a  generalized  phase  space 
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<>> 


where  all  time  derivatives,  save  for  the  one  of  highest  order,  become  independent  variables  [3]. 
A  generic  linear  homogeneous  first-order  system  can  be  written  in  the  form 

where  a  linear  operator  'H  is  called  Hamiltonian,  while  ’J't  is  called  a  state  vector  (or  wave 
function)  .  It  is  an  element  of  the  generalized  phase  space  of  the  system  and  viewed  a  collection 
(column)  of  the  original  variables  and  their  time  derivatives.  The  generalized  phase  space  is 
equipped  with  an  inner  product  and  becomes  a  Hilbert  space.  State  vectors  are  typically 
vector-valued  functions,  and  the  Hamiltonian  is  a  differential  operator.  The  choice  of  the 
inner  product  depends  on  the  problem  on  hands.  One  usually  requires  that  components  of 
are  elements  of  the  space  of  square  integrable  functions. 

The  advantage  of  the  Hamiltonian  formalism  is  evident  because  the  fundamental  solution 
can  be  written  as  the  exponential  of  the  Hamiltonian, 

=  exp  im)  ’J'o  =  >  0  >  (1*^) 

assuming,  of  course,  that  the  exponential  of  U  exists.  In  numerical  studies,  the  exponential 
representation  of  the  fundamental  solution  yields  a  rather  natural  time  domain  algorithm. 
The  amplification  matrix  Wa*  can  effectively  be  computed  by  means  of  the 
Kato-Trotter  product  formula  in  combinations  with  pseudospectral  methods  [2]. 

The  Hamiltonian  formalism  has  been  used  in  [4]  to  develop  a  finite  differencing  algorithm 
to  study  an  electromagnetic  pulse  propagation  in  Lorentz  media.  Consider  2N  real  vector 
fields,  ^1,  j  =  1, 2, ...,  2N,  such  that 

p;  =  ,  (1-7) 

Thus,  the  original  system  of  second  order  equations  has  been  converted  into  the  first  order 
system.  Another  convenient  way  to  introduce  the  Hamiltonian  formalism  is  to  use  N  complex 
vector  fields  Ct  which  satisfy  the  first  order  differential  equation 

dtCt  =  ,  P?  =  ^  (Ct“  +  Ct )  ,  (1-9) 

where  Aa  =  -7a  +  and  i/*  =  y'u;2  -  72.  This  representation  is  defined  only  if  7a  <  Wo  (i.e., 
the  attenuation  is  not  high).  From  the  numerical  point  of  view,  solving  a  decoupled  system  of 
N  first  order  differential  equation  and  taking  complex  conjugation  (denoted  here  by  an  over 
bar)  is  less  expensive  than  solving  an  original  system  of  differential  equations  to  compute  the 


*'=(^<),  (1.10) 


medium  polarization. 
Define 


•■'-(t')- 
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where  the  vector  is  composed  of  and 

<=©’  '^•  =  (b:)' «=(tr)- 

The  wave  function  satisfies  the  Schrodinger  equation 


where 


(  ‘Ho 
VVmf 

Vfm 

'tjF 

%0 

/  0  cVx\ 

\-cVx  0  ) 

5 

Vfm 

= 

(VfMIj  VfM2i  • '  •  >  Vfmn)  ) 

VfMo  ■ 

II 

Vmf 

= 

-Vfm  , 

-t/F 

diag  {'Hmij  'Hm2i’ 

'  *  5  H-mn) 

,  H 

F  _  1 
Ma  —  1 

’ 


0  )  ’ 


0 

-Wa 


(^a 

-H 


fa)  ’ 


(1.11) 


(1.12) 

(1.13) 

(1.14) 

(1.15) 


where  diag  indicates  that  the  corresponding  matrix  i&  block-diagonal  with  blocks  listed  in  the 
order  from  the  upper  left  to  lower  right  comers.  Note  that  the  matrices  VpMa  and  HpMa  act 
on  a  six-dimensional  column  Therefore  they  should  be  understood  as  composed  of  3  x  3 
blocks.  Each  block  is  obtained  by  multiplying  the  unit  matrix  by  the  number  indicated  in 
place  of  the  block  in  (1.13)  and  (1.15).  We  shall  refer  to  (1.11)  as  to  the  theory  in  the  field 
representation. 

If  the  inductions  are  used  as  independent  electromagnetic  degrees  of  freedom,  this  is, 
the  corresponding  evolution  equation  can  be  obtained  by  the  similarity  transformation.  Let 

be  an  operator  that  maps  the  auxiliary  matter  field  onto  the  actual  medium  polarization 
vector,  Pt  =  Its  explicit  form  is  not  hard  to  deduce.  Then 


.;=5«i»r,  5=(iT)- 

sn^s-^ . 


f)’ 


Hence, 

The  corresponding  blocks  of  have  the  form 


‘Hi  —  "Ho  ,  Vmi  =  VftfF  5 
Vim  =  Vfm  +  'H'Hm  “■  'Ho'R  —  —'Ho'R  > 

njl  _ 

TtnA  — 


M-M  —  ■ 


(1.16) 

(1.17) 

(1.18) 

(1.19) 

(1.20) 


In  what  follows  we  shall  use  where  Q  indicates  the  representation,  F  or  7. 
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Our  final  remark  concerns  a  change  of  basis  in  the  space  of  auxiliary  fields  (canonical 
transformations  in  the  Hamiltonian  formalism).  There  is  a  freedom  of  making  general  complex 
nonsingular  linear  transformations 

det^M^O.  (1.21) 


If  the  infinitesimal  evolution  operator  U^t  =  exp(AtK^)  is  computed  with  one  choice  of  the 
auxiliary  fields,  a  simple  similarity  transformation,  like  the  one  in  (1.17),  would  allow  us  to 
compute  it  in  any  other  basis  of  auxiliary  fields.  This  is  an  important  observation  because 
the  auxiliary  field  basis  can  be  chosen  in  a  way  that  facilitates  computation  of  the  evolution 
operator  (e.g.,  to  improve  the  convergence  rate  or  speed  up  simulations).  For  instance,  in 
the  complex  representation  (1.9)  of  the  auxiliary  fields  in  the  Lorentz  model,  the  matter 
Hamiltonian  is  diagonal.  The  corresponding  transformation  of  the  auxiliary  fields  is  given 


To  transform  the  whole  system  into  this  representation,  the  Hamiltonian  is  replaced  by 
and  the  wave  function  'J'f  by  where  S  is  block-diagonal  with  the  unit  matrix 
in  the  upper  left  (field)  comer  and  with  Sm  in  the  lower  right  (matter)  comer. 


2  Energy  and  the  norm  of  state  vectors 


Let  us  discuss  the  choice  of  a  norm  in  the  space  spanned  by  state  vectors  Consider 
multi-resonant  Lorentz  models  with  no  attenuation  'ya  =  0.  The  field  and  matter  evolution 
equations  can  be  obtained  from  the  variational  principle  for  the  action 


=  J  dtL  =  I dtj dr  i  (B?  -  B?)  +  i  ^ 


.  (2-1) 


where  electromagnetic  degrees  of  freedom  are  described  by  the  vector  and  scalar  potentials, 
respectively.  At  and  (ft,  so  that  Et  =  VyJt  ~  dtA.t  and  Bt  =  V  x  At.  The  units  are  chosen 
in  this  Section  so  that  c  =  1.  The  polarization  of  the  medium  is  expressed  via  the  matter 
fields  as  Pt  =  action  principle  for  the  scalar  potential  (pt  leads  to  the 

Gauss  law,' V  -Dt  =  0,  for  the  vector  potential  At  to  the  Maxwell’s  equation,  ^Dt  =  V  x  Bt, 
and  for  the  matter  field  t?®  to  the  medium  polarization  evolution  equation  of  the  Lorentz 
model  with  no  attenuation,  ja  —  0*  The  second  Maxwell’s  equation  and  the  Gauss  law  for 
the  magnetic  field  follows  from  the  relation  Bt  =  V  x  At  by  taking  its  time  derivative  and 
divergence,  respectively.  The  energy  of  the  system  coincides  with  the  canonical  Hamiltonian 
which  is  obtained  by  the  Legendre  transformation  [3]  of  the  Lagrangian  L  for  the  velocities 
dtAt  and  dtOf.  Doing  the  Legendre  transform,  we  find  the  canonical  Hamiltonian  (energy)  of 


the  system 


+  B^,+-£{nf  +  u,l0f) 

a 


(2.2) 


5 


». 


where  tt^  =  5L/5{dt‘d°')  are  canonical  momenta  of  the  matter  fields.  The  energy  conservation 
follows  directly  from  the  Noether  theorem  [3]  applied  to  the  time  translation  symmetry  of  the 
action  (2.1),  =  0.  Equation  (2.2)  becomes  the  conventional  expression  for  the  electro¬ 

magnetic  energy  in  a  passive  medium  [1]  when  7r“  and  are  replaced  by  the  corresponding 
solutions  of  the  equations  of  motion  with  initial  conditions  Trg  =  =  0. 

An  important  observation  is  that  the  Noether  integral  of  motion  (2.2)  coincides  with  the 
norm  squared  of  the  corresponding  state  vector  [2] 

Et  =  ^J  dr  =  ('i^f ,  )  =  (^^,  ,  (2.3) 

where  fi  =  This  is  proved  by  observing  that  which 

follow  from  comparison  of  the  canonical  Hamiltonian  equations  of  motion  [3]  for  the  canon¬ 
ically  conjugated  variables  ^  and  tt®  and  Eqs.  (1.8)  with  7®  =  0.  Note  that  the  canonical 
momentum  conjugated  to  the  vector  potential  Af  coincides  Avith  — Dt  =  — E*  —  Pt,  not  — Ef 
in  this  system.  Therefore,  the  coupling  between  the  electromagnetic  and  matter  degrees  of 
freedom  is  included  into  the  term  E?  =  (Dt  -  Pt)^  of  the  canonical  Hamiltonian  (2.2). 

Thus,  in  the  absence  of  attenuation,  the  norm  of  the  state  vector  is  proportional  to  the 
wave  packet  electromagnetic  energy  and,  hence,  is  conserved.  The  norm  conservation  follows 
also  from  the  skew-symmetry  of  the  Hamiltonian  'H^*  =  if  7o  =  0)  while  (2.3)  establishes 
a  relation  between  the  electromagnetic  energy  and  the  norm.  In  the  induction  representation, 
the  norm  in  the  measure  space,  defined  by  the  operator  fi  in  (2.3),  is  also  conserved  by 
construction.  Consequently,  the  Hamiltonian  is  skew-symmetric  relative  to  the  measure  space 
scalar  product, 

The  norm  (energy)  conservation  can  be  used  to  control  numerical  convergence,  especially 
when  the  aliasing  problem  in  the  fast  Fourier  transform  is  present,  e.g.,  when  parameters 
of  the  medium  are  discontinuous  functions  in  space.  In  a  properly  designed  algorithm  the 
loss  of  energy  (norm)  due  to  attenuation  should  be  controlled  by  the  symmetric  part  of  the 
Hamiltonian  operator 

=  V«*?)<0,  (2.4) 

a 

where  V^*  =  =  {'HP*  +  H*^)/2  <  0  (a  negative  semidefinite  operator)  which  is,  in  this 

case,  a  diagonal  matrix  with  nonpositive  elements. 

If  the  medium  in  question  does  not  have  dispersion  and  absorption,  the  formalism  is  sim¬ 
plified.  Let  e  =  e(x)  be  the  dielectric  constant  of  the  medium.  If  the  medium  is  not  isotropic, 
then  e  is  symmetric  positive  definite  3x3  matrix  everywhere  in  space.  We  rewrite  Maxwell’s 
equations  in  the  form 

dtipt  =  ^  j  0  )  ’ 
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where  the  brackets  in  (e~^)  mean  that  the  induction  is  first  multiplied  by  e  and  then  the 
curl  of  the  resulting  vector  field  is  computed.  Consider  the  scalar  product 

=  (^0  l)' 

The  Hamiltonian  is  anti-Hermitian  with  respect  to  this  scalar  product,  MqHs  = 
Therefore  the  corresponding  fjt  norm  is  preserved  in  the  time  evolution  generated  by  ejqpitHa), 
that  is,  The  electromagnetic  energy  of  the  wave  packet  is  conserved 

because  it  is  proportional  to  the  fj,  norm  of  the  initial  state  vector. 


3  The  algorithm 

Consider  an  equidistantly  spaced  finite  grid  with  periodic  boundary  conditions.  Let  Ar  be  the 
grid  step  and  n  be  a  vector  with  integer  valued  components.  Then  the  dynamical  variables 
are  projected  onto  the  grid  by  taking  their  values  at  grid  points  r  =  iiAr,  that  is,  (r)  -> 
^f{nAr).  A  cubic  grid  is  assumed.  Consider  a  discrete  Fourier  transformation  associated 
with  the  grid,  W?(nA:o)  =  X)n' where  =  .FP  =  1.  The  dual  lattice 
step  is  ko  =  27r/Ar.  The  grid  spatial  size  L  and  step  Ar  must  be  chosen  so  that  the  Fourier 
transform  of  the  initial  wavepacket  has  support  within  the  region  k  G  where 

ik  =  |k|,  kmax  =  h  and  kmin  =  2tc/L.  The  Hamiltonian  is  split  into  a  sum  'HP  =  'Hq  +  V^, 
where  all  the  spatial  derivatives  are  included  into  ‘H!q  and  contains  multiplications  by 
position  dependent  functions.  The  operator  is  projected  naturally 

(r)  V^(nAr)^f  (nAr)  .  (3.1) 

Consider  Hq  in  the  Fourier  basis,  'Ho{ik).  The  projection  is  then  done  via  the 

discrete  Fourier  transform 

H?(V)»?(r)|  .  (3.2) 

I  r=nAr 

n' 

The  projection  (3.2)  as  well  as  any  action  of  'Hq  on  state  vectors  is  performed  by  the  fast 
Fourier  method. 

In  what  follows,  the  rules  (3.1)  and  (3.2)  define  the  action  of  the  operators  and  'Hq 
and  their  functions  on  any  state  vector.  The  action  of  a  product  of  and  Hq  on  any  state 
vector  is  understood  as  consecutive  actions  of  these  operators  according  to  the  rules  (3.1)  and 
(3.2)  in  the  order  specified  in  the  product.  The  main  advantage  of  using  the  Fourier  basis  is 
the  exponential  convergence  (versus  the  polynomial  one  in  finite  differencing  schemes)  [5]  as 
the  grid  size  increases. 

The  time  evolution  on  the  grid  is  simulated  by  the  time  leapfrog  method  for  the  Schrodinger 
equation  (1.5)  in  both  the  cases,  the  Lorentz  model  and  the  geometric  optic  case, 

%+At  =  %-At  +  2AtH%  .  (3.3) 
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As  is  well  known  from  the  Fourier  analysis,  the  convergence  rate  can  be  affected  for  functions 
which  have  discontinuities  [6].  The  latter  is,  unfortunately,  the  case  in  electromagnetic  scat¬ 
tering  problems.  Suppose  there  is  an  interface  between  two  media.  It  can  be  deduced  from 
the  Maxwell’s  equations  that  the  components  of  the  electric  and  magnetic  fields,  Et  and  H*, 
tangential  to  the  interface  must  be  continuous,  provided  there  is  no  surface  electric  current 
on  the  interface.  From  the  Gauss  law  it  follows  that  the  components  of  the  inductions,  Dt 
and  Bt,  normal  to  the  interface  must  be  continuous,  provided  there  is  no  surface  charge  on 
the  interface.  In  contrast,  the  normal  components  of  the  fields  and  the  tangential  components 
of  the  inductions  can  be  discontinuous.  Their  discontinuities  are  proportional  to  discontinu¬ 
ities  of  medium  parameters  (e.g.,  discontinuities  in  plasma  frequencies  in  Lorentz  models). 
Therefore,  in  either  the  induction  or  field  representation,  there  are  components  which  suffer 
discontinuities  at  the  interface.  The  only  way  to  cope  with  the  problem,  while  keeping  the 
use  of  the  Fourier  basis,  is  to  make  the  grid  finer.  This  would  lead  to  a  substantial  waste 
of  computational  resources  because  the  sampling  efficiency  should  only  be  increased  in  the 
neighborhood  of  medium  interfaces.  The  use  of  wavelet  bases  might  be  efficient  for  such  a 
task  [7].  Here  we  retain  the  Fourier  basis,  and  increase  the  sampling  efficiency  by  a  change  of 
variables. 

Consider  a  change  of  variables  y  =  y(x).  A  uniform  grid  in  the  new  variables  y  would 
generate  a  non-uniform  grid  in  the  original  Euclidean  (physical)  coordinates  x.  A  desired  local 
density  of  grid  points  in  the  physical  space,  to  enhance  the  sampling  efficiency  in  designated 
regions,  can  be  achieved  by  an  appropriate  choice  of  the  functions  y(x)  [8].  The  fast  Fourier 
transform  algorithm  is  applied  to  the  auxiliary  uniform  grid  according  to  the  rule  V*  = 
A{y)Vy,  where  the  product  of  the  derivative  and  multiplication  operators  is  treated  by  the 
rules  (3.2)  and  (3.1). 

Let  us  introduce  the  amplification  matrix  for  the  leapfrog  algorithm,  ’i't+At  =  It 

satisfies  the  equation  Git  —  2AtHGAt  —  1  =  0  which  has  two  solutions 

=  HAt  ±Vl+  .  (3.4) 

The  stability  of  the  algorithm  requires  that  the  energy  norm  of  both  the  solutions  (^S^)”’®^o 
must  be  bounded.  The  necessary  condition  (but  not  sufficient)  is  the  von  Neumann  condition 
that  the  spectral  radius  of  the  amplification  matrix  does  not  exceed  1.  If  is  an  eigenvalue 
of  AtH,  then  the  von  Neumann  condition  leads  to  ip  =  ±7r/2  and  <  1.  In  other  words,  the 
eigenvalues  must  be  imaginary  and  their  magnitude  should  not  exceed  1.  So,  the  Hamiltonian 
should,  at  least,  be  anti-Hermitian,  and  therefore  the  von  Neumann  criterion  becomes  also 
sufficient  for  the  stability  because  anti-Hermitian  operators  are  also  normal. 

In  the  geometric  optic  case,  the  Hamiltonian  is  anti-Hermitian,  and  the  von  Neumann 
condition  is  fulfilled  if 

1  > 

where  kt,„^  the  maximal  norm  of  all  wave  vectors  in  the  medium  which  can  be  estimated  by 
\/p{^)kmax  with  kmax  being  the  maximal  wave  vector  of  the  initial  pulse  in  vacuum  and  p{e} 
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the  maximal  spectral  radius  of  the  symmetric  matrix  e{x)  over  x  (or  simply  the  maximum 
of  e{x)  if  the  medium  is  isotropic).  This  can  be  understood  from  the  following  pnnciple  [9J. 
A  finite  difference  scheme  with  variable  coeflacients  is  stable  if  all  the  corresponding  schemes 
with  frozen  (i.e.,  fixed  to  a  particular  value  everywhere  in  space)  coefficients  are  stable. 


The  von  Neumann  condition  cannot  be  met  for  Lorentz  models  if  absorption  is  present, 
7a  ^  0.  To  circumvent  this  difficulty,  the  leapfrog  scheme  is  modified  in  the  following  way 
[2]  Letn  =  no  +  V  where  Uq  =  -'Hq  is  such  that  the  von  Neumann  condition  is  satisfied 
for  it,  while  V  is  negative  semidefinite,  that  is,  for  any  Vt)  £  0.  In  Lorentz  models 

this  is  easily  achieved  by  taking,  for  instance,  Mo  =  'H\^a=o-  Then  V  is  diagonal  with  matrix 
elements  being  zeros  and  -27a.  Another  possibility  is  to  identify  Uq  with  the  Hamiltonian  in 
the  vacuum,  then  it  is  straightforward  to  show  that  V  =  %- Ho  is  negative  semidefinite.  In 
(1.5)  we  make  a  substitution  =  exp(tV)$t.  The  new  state  vector  satisfies  the  equation 
with  a  time  dependent  Hamiltonian 

=  e-*^noe*'^^t  =  ,  (3-5) 


and  with  the  same  initial  condition  $o  =  Applying  the  leapfrog  method  to  (3.5)  we  get 
^t+Af  =  ^t-At  +  valid  up  to  0(At®).  Returning  to  the  initial  variables,  we  arrive  at 

the  following  recurrence  relation 

'i't+At  =  ,  (3-6) 


where  Cm  =  exp(AtV). 
satisfies  the  equation 


The  amplification  matrix,  ’J^t+At  —  for  recurrence  (3.6) 

Qm  —  C2mQm  2AtCMHo  .  (3*7) 


According  to  Theorem  7  of  [2]  the  deviation  of  the  approximate  solution  fro“i  the  exact 

solution  relative  to  the  energy  norm  is  of  order  At^  for  any  n.  Thus,  a  possible  exponential 
norm  growth  can  be  suppressed  as  much  as  desired  by  making  the  time  step  small  enough. 
Mathematical  details  of  the  proof  are  omitted  here.  However,  its  basic  idea  can  be  understood 
from  the  following  observation.  Solving  (3.7)  by  perturbation  theory  in  At,  it  is  not  hard  to 

find  that 

Qm  ~  Qm  ~  A^K.m  >  Q\t  ~  CmiiQ •>  ^ 


where  Km  is  regular  in  the  vicinity  of  At  =  0  and  vanishes  whenever  Uq  and  V  commute, 
and  Q^t  is  amplification  matrix  when  V  is  set  to  zero.  The  von  Neumann  condition  is 
satisfied  for  Qm  (or,  it  is  always  possible  to  find  At  >  0  to  satisfy  it).  Hence  the  lowers  of 
Q%t  applied  to  ’J'o  do  not  produce  the  exponential  norm  growth.  The  powers  of  ^Xt  differs 
from  those  of  by  factors  which  can  only  produce  the  exponential  attenuation  of  the 
amplitude.  Indeed,  let  Then  =  2(’3?r>V^n  <  0  since  V  is  negative 

semidefinite.  Thus,  the  solution  produced  by  the  amplification  matrix  ^Xt  no  exponential 
growth,  while  differing,  relative  to  the  energy  norm,  from  that  produced  by  Q^t  by  order  of 
O(Ai^).  Since  both  the  solutions  differ  from  the  exact  one  by  order  of  0{Af),  we  conclude 
that  the  modified  leapfrog  scheme  can  be  made  as  accurate  as  desired  by  reducing  the  time 
step.  Here  is  an  example. 
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In  the  field  representation  of  the  Hamiltonian  for  multiresonant  Lorentz  models,  we  make 


the  following  decomposition 

V)  +  (o 

As  a  result  we  arrive  at  the  following  scheme 

=  V’f-At  +  2At?{o^'^  +  2AiX;>^FMa?t  ,  (3-10) 

a 

C+At  =  +  2At  VMPa  i’t  >  {3-11) 

e*«Mo  =  g-T'o*  cosh  i>at  +  —  {'Hmu  +  7o)  >  (3.12) 


where  i>a  =  The  exponential  (3.12)  is  easy  to  compute  by  expanding  in 

the  Pauli  matrix  basis,  which  is  also  a  basis  for  the  Lie  algebra  su(2),  and  then  by  using 
the  well  known  formula  for  the  exponential  of  a  linear  combination  of  Pauli  matrices.  For 
small  attenuation,  7„  <  Wa,  we  get  i>a  =  iva-  The  hyperbolic  functions  in  (3.12)  become  the 
trigonometric  ones  and  Va  is  replaced  by  The  eigenvalues  of  the  matter  Hamiltonian  are 
^  ±  Hence,  Re  <  0  and  amplitudes  of  the  matter  fields  are  always  exponentially 

attenuated  as  t  oo,  unless  7®  =  0  leading  to  Re  Ao  =  0. 

The  stability  is  ensured  if  satisfies  the  von  Neumann  condition.  Let  kmax  be  the 
maximal  norm  of  all  wave  vectors  of  the  initial  wave  packet  and  be  the  maximal  value 
of  U}p  =  as  a  function  of  position,  then  a  sufficient  criterion  for  stability  reads 

(3-13) 

Here  the  idea  of  the  frozen  coefficients  has  been  used  again.  The  scheme  (3.6)  becomes 
especially  simple  in  the  case  of  small  attenuation  7a  <  In  the  complex  representation  of 
the  auxiliary  fields  (1.22)  (cf.  (1.9))  the  matter  Hamiltonians  are  diagonal  and  the  action 
of  its  exponential  is  reduced  to  multiplication  by  a  complex  number 

Finally,  it  should  be  mentioned  that,  by  rearranging  operators  in  the  split,  namely,  by 
moving  Vfm  and  Vmf  to  in  (3.9),  the  stability  condition  (3.13)  can  be  weakened  to 
Atckmax  <  1*  This  would  come  at  the  price  of  having  a  more  complicated  expression  for  C^t- 
In  the  case  of  the  Lorentz  model  it  can  still  be  computed  analytically.  The  new  split  can 
also  be  viewed  as  the  use  of  the  induction  representation  in  the  modified  leapfrog  scheme, 

= 'Hi where  'Hi  contains  only  the  blocks  of  'H^  with  the  derivative  operator  V.  The 
proof  of  the  weaker  stability  condition  can  be  found  in  [2]. 
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